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Summary 


A study has been performed focusing on the calculation of sensitivities of displacements, 
velocities, accelerations, and stresses in linear, structural, transient response problems. One 
significant goal of the study was to develop and evaluate sensitivity calculation techniques 
suitable for large-order finite element analyses. Accordingly, approximation vectors such as 
vibration mode shapes are used to reduce the dimensionality of the finite element model. Much 
of the research focused on the convergence of both response quantities and sensitivities as a 
function of the number of vectors used. 

Two types of sensitivity calculation techniques were developed and evaluated. The first type 
of technique is an overall finite difference method where the analysis is repeated for perturbed 
designs. The second type of technique is termed semianalytical because it involves direct 
analytical differentiation of the equations of motion with finite difference approximation of the 
coefficient matrices. To be computationally practical in large-order problems, the overall finite 
difference methods must use the approximation vectors from the original design in the analyses 
of the perturbed models. This was found to result in poor convergence of stress sensitivities 
in several cases. To overcome this poor convergence, two semianalytical techniques were 
developed. The first technique accounts for the change in eigenvectors through approximate 
eigenvector derivatives. The second technique applies the mode acceleration method of transient 
analysis to the sensitivity calculations. Both result in very good convergence of the stress 
sensitivities. In both techniques the computational cost is much less than would result if the 
vibration modes were recalculated and then used in an overall finite difference method. 



Chapter 1 
Introduction 


1.1. Overview 

In the past 10 years there has been increasing in- 
terest in calculating the derivatives of structural be- 
havior with respect to problem parameters or design 
variables (i.e., sensitivities). One of the main uses of 
these sensitivities is in automated design procedures 
where a numerical algorithm is used to improve a 
structure by modifying the design parameters while 
satisfying prescribed constraints on the structural be- 
havior. Most of the numerical algorithms used in 
these procedures require both an initial design and 
a set of sensitivities in order to decide how to im- 
prove the structure. Many references address this 
sensitivity calculation question within the context of 
automated structural design, whereas others, such as 
this study, focus specifically on issues related to the 
calculation of sensitivities. Other uses of sensitivities 
in structures problems include the system identifica- 
tion problem in structural dynamics and statistical 
structural analysis. References 1 and 2 provide a 
comprehensive review of work on calculating sensi- 
tivities in structural systems. 

It is clear from many references (e.g., ref. 1) 
that most of the emphasis in structural optimiza- 
tion and the associated sensitivity calculation meth- 
ods has been on static problems. This is not sur- 
prising since most structural analyses themselves are 
static. The objective of the static analysis and sen- 
sitivity calculation problem, for linear systems, is to 
calculate the responses (e.g., displacements, stresses) 
and their derivatives with respect to structural pa- 
rameters (e.g., member areas, thicknesses), which 
are assumed to be constant for all time. Tech- 
niques for both the analysis and sensitivity calcula- 
tions have reached considerable maturity in the past 
10 to 20 years. 

In many problems, however, the loading on the 
structure varies with time, which causes the response 
of the structure also to vary as a function of time. 
Examples of such problems are a gust on an aircraft 
wing, an unbalanced engine in an automobile, or a 
building during an earthquake. In these cases, it is 
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important to predict stresses accurately as well as 
displacements (and possibly velocities and accelera- 
tions) as a function of time. Often it is sufficient to 
predict the maximum and minimum values of these 
response quantities. Similarly, the goal of the sen- 
sitivity analysis is the calculation of derivatives of 
these response quantities with respect to the struc- 
tural parameters as a function of time or at the time 
points where the maximum or minimum responses 
occur. 

The introduction of the time parameter compli- 
cates the analysis in several ways. First, it changes 
the system of equations from a set of coupled alge- 
braic equations to a set of coupled differential equa- 
tions whose accurate solution may be difficult and 
computationally costly. Second, the amount of in- 
formation that must be considered and evaluated to 
understand the response of the structure is increased 
by orders of magnitude. 

Most practical static and dynamic analyses are 
currently performed with the finite element method. 
Since this technique replaces a continuum (infinite di- 
mensional space) with a finite-degree-of- freedom ap- 
proximation, the question of required mesh refine- 
ment is a natural one. This is not an easy question to 
answer because the convergence of the approximation 
as the mesh is refined depends on the quantity being 
considered. Usually, the fundamental unknowns are 
the displacements and rotations at the finite element 
nodes. In theses cases, the convergence of derivatives 
of displacements with respect to a spatial parameter 
(stresses), with respect to time (velocities, acceler- 
ations), or with respect to a structural parameter 
(sensitivities) will be worse than the convergence of 
the displacements themselves. 

After the structure has been discretized with the 
finite element method, yet another approximation 
is usually introduced in linear dynamics problems. 
The behavior of the structure is represented by a 
reduced set of basis functions (frequently natural 
vibration modes) in order to simplify the solution of 
the transient response problem. This approximation 
introduces another set of concerns over accuracy of 
the response quantities and their sensitivities. 

Other errors in transient analysis or sensitivity 
calculations, which rarely occur in static analyses, 
are due to the truncation error of finite difference 
operators. This problem occurs with the use of nu- 
merical integration techniques in solving the coupled 
differential equations in the transient problem. This 
problem also occurs when difference approximations 
are used in the calculation of sensitivities. Round-off 
errors, due to the finite precision arithmetic on digi- 
tal computers, are also more of a concern in transient 
or sensitivity analyses than in simple static analyses. 
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All these complexities in transient analysis cou- 
pled with the problems of sensitivity analysis have 
slowed the progress in the development of sensi- 
tivity calculation techniques for transient response 
problems. However, substantial progress has been 
made. Some of the important previous work in 
optimization of structures under transient loads and 
calculation of sensitivities in transient response prob- 
lems is discussed in section 1.2. 

1.2. Review of Previous Pertinent Work 

Reference 3 is one of the earliest papers dealing 
with optimization of structures under transient loads. 
In this paper, Fox and Kapoor consider the minimum 
mass design of frame structures under an applied 
base motion subject to constraints on deflections and 
stresses. The equations of motion are uncoupled by 
using vibration modes and solved for the maximum 
value of the modal response by using a shock spec- 
trum approach. A considerable simplification is in- 
troduced by directly summing the maximum modal 
responses; therefore, time is removed as a parameter 
in the calculations. 

In references 4 and 5, Cassis and Schmit present 
procedures for the automated design of plane frames 
under general transient loading. The dynamic analy- 
sis is performed with modal superposition, and only 
modal damping is allowed. Integral forms of the 
time-dependent constraints are used. Sensitivities 
are calculated with an explicit differentiation of the 
dynamic equations along with exact calculation of 
the required eigenvalue and eigenvector derivatives. 
Effects of finite element discretization and modal 
truncation on the sensitivities or final optimized de- 
signs were not considered. 

In the past 10 years, other researchers have con- 
sidered the application of general sensitivity theory 
to the problem of dynamic mechanical systems. Ref- 
erence 1 summarizes this work and describes three 
basic approaches which have been employed. In the 
first method, called the direct method, the equations 
of motion are directly differentiated and solved. A 
second method offers the advantage of reduced com- 
putational cost when there are more design vari- 
ables than constraints on response quantities. In 
this method, called the adjoint method, the sensi- 
tivity equations are rewritten in terms of a newly 
defined adjoint vector. After solving this new sys- 
tem for the adjoint vector, the calculation of the sen- 
sitivities of the response constraints with respect to 
each of the design variables is straightforward. In the 
third method, called the Green’s function method, 
the derivatives are obtained in terms of the Green’s 
function of the equations of motion. Although the 


results from all three methods are theoretically iden- 
tical, their relative computational efficiency depends 
on the relative numbers of design variables, degrees 
of freedom, and constraints. 

Haug, Arora, and their coworkers have made con- 
siderable progress in addressing many of the prob- 
lems in the optimal design of mechanical systems un- 
der dynamic loadings. Much of their early work was 
spent studying a “state space” or adjoint variable 
approach to calculating sensitivities. References 6 
through 9 should be noted. These references con- 
sider application to both elastic structural design and 
machine design problems that often have the addi- 
tional complexity of nonlinear equations of motion. 
However, most of these examples have involved few 
degrees of freedom or design variables. A more re- 
cent paper by Haug (ref. 10) extended the sensitivity 
analyses of previous papers to include additional al- 
gebraic constraint equations that are often present in 
machine design problems. Also, sensitivity equations 
for second derivatives are presented. 

The adjoint method is particularly attractive 
when a transient constraint is integrated over time 
to produce a single constraint because the total num- 
ber of constraints is often small. However, the loss 
of information in this integral formulation and its 
disadvantages are noted in reference 11. Given the 
danger of having only a single “worst-case” value of 
the constraint function in time, reference 11 proposed 
including all local maximum points of the constraint 
function in the constraint set. A significant disadvan- 
tage of this approach is that for “jagged” response 
functions, there can be a large number of redundant 
local maxima. This important problem of constraint 
definition was also considered in reference 8, where 
several methods for obtaining a few important con- 
straints at discrete points in time were proposed. 

Both direct and adjoint sensitivity methods for 
a nonlinear hysteretic structure are presented in ref- 
erence 12. Because of the nonlinearities, numerical 
integration of the full coupled system is required. 

A recent approach in sensitivity analysis has been 
to write sensitivity expressions for the solid contin- 
uum prior to discretizing the system. This approach 
is especially attractive when shape-type design vari- 
ables are being considered because the design vari- 
able itself often represents a continuous region on 
the surface of the body. Reference 13 uses the con- 
cept of the material derivative to calculate shape 
derivatives of a continuum under dynamic loads. In 
reference 14, expressions for shape sensitivities of 
a continuum considering material nonlinearities and 
dynamic effects are written with a variational 
approach. 


4 



1.3. Objectives and Scope 

The purpose of the study reported herein is to in- 
vestigate methods for calculating sensitivities in lin- 
ear transient structural response problems. Very gen- 
eral forms of external loading on the structure and 
damping are permitted. In any numerical algorithm, 
both accuracy and computational efficiency are con- 
cerns. Errors in the sensitivities due to factors such 
as the finite element mesh, truncation of the basis 
vector set in the transient analysis, and finite differ- 
ence approximations in the sensitivity and numerical 
integration procedures are considered. An objective 
of the study is to identify approaches to sensitivity 
analysis that are appropriate for large-scale struc- 
tural analysis. This is emphasized in the selection of 
the algorithms and in a study of the relative compu- 
tational efficiency of several competing methods. 

Three transient response problems are considered 
in detail: a five-span, simply supported beam; a com- 
posite aircraft wing; and a cantilever beam with a 
cross section that varies along its length. None of 
these three problems are large. However, each prob- 
lem includes ingredients which make the sensitivity 
analysis computationally difficult. 




Chapter 2 

Equations of Motion and Solution 

2.1. Governing Equations 

The equations of motion for a damped, linear 
structural system can be written as 

Mii + Cu + Ku = p (t) (2.1) 

which is a set of n g coupled differential equations 
and M, C, and K are the system mass, damping, 
and stiffness matrices, respectively. Frequently it 
is possible to separate the loading vector p into a 
product of a vector describing the spatial distribution 
of the loading f and a scalar function of time g{t) as 

P (t) = g(t ) f (2.2) 

Often equations (2.1) are the result of a large 
finite element model and are therefore of large order. 
One way to characterize the behavior of this system 
is by examining the eigenvalues of the undamped 
system 

K</»j - wJm 4>j = 0, = rig) (2.3) 

For most large structural systems, equations (2.1) are 
“stiff” ; the condition number /u is many orders 
of magnitude. 

The external loading also has a major effect on 
the dynamic response of the system. Impulsive loads 
where g(t) changes rapidly relative to the periods as- 
sociated with the smallest ujj tend to produce a re- 
sponse history with significant high-frequency com- 
ponents. Loads that are applied slowly relative to the 
vibration periods of the c oj produce a predominantly 
low-frequency response history. 

Two basic approaches are available for the solu- 
tion of equations (2.1). The first approach is to nu- 
merically integrate the equations in a step-by-step 
manner. In implicit integration techniques, the time 
step must be a fraction of the period associated with 
the largest t oj significantly excited by the loading in 
order to obtain an accurate solution. The well-known 
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Ncwmark method (for example, ref. 15) is an example 
of such an integration technique. In explicit integra- 
tion techniques, the time step must be a fraction of 
the period associated with u) Ug in order for the solu- 
tion process to be numerically stable. Using either 
technique, the computational work is large because 
equations (2.1) are of large order. 

An alternative to directly solving equations (2.1) 
is to solve an approximate reduced-order problem 
instead. This is the preferred approach for most 
linear structural dynamics problems. The details 
of the techniques used to reduce the order of the 
dynamic system are discussed in section 2.2. 

2.2. Reduction Techniques 

The first step in applying a reduction technique 
to the solution of equations (2.1) is to approximate 


the solution by n r basis functions 

u = 4>q (2.4) 

where n r is usually much less than n g . Then a 
reduced set of equations can be written 

Mq A Cq + Kq = g{t)f (2.5) 

where 

M = (2.6) 

C = # T C$ (2.7) 

K = <f> T K$ (2.8) 

f = <fc r f (2.9) 


If the number of vectors in is equal to the size of the 
original system n g and the vectors in are linearly 
independent, the transformation of equation (2.4) is 
exact. Usually, though, n r « n g and the solution 
to the full system (eqs. (2.1)) is only approximated 
by the solution to the reduced system (eqs. (2.5)). 
The quality of this approximation as the number 
of vectors in 4> is increased is a key concern in 
evaluating the effectiveness of a particular reduction 
technique. 

In all reduction methods considered herein, the 
first n r vectors of the set are taken as the reduced 
basis. Alternate approaches are available for assess- 
ing the importance of a given vector prior to solution 
of the reduced system and then discarding the vector 
if its contribution is insignificant. These approaches 
are not considered here because the cost of generat- 
ing the set of vectors is often high and the cost of 
solving equations (2.5) is often fairly low. 
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2.2.1. Mode Displacement Method 

The most widely used reduction technique is 
the traditional mode displacement method. In this 
method, equations (2.3) are solved for the set of vi- 
bration modes with lowest n r frequencies and modes. 
This set of vibration modes is used as the set of basis 
functions When the system is undamped (C = 0 
in eqs. (2.1)) or C can be expressed as a linear combi- 
nation of M and K, equations (2.5) represent a set of 
uncoupled differential equations which can be solved 
independently. If the eigenvectors are scaled so that 
(f)jNL(f)j — 1, the uncoupled equations can be written 
as 


( ii + 2 ijUiqi + Jjm - y(t)f, (i = 1 , — , n r ) 

( 2 . 10 ) 

where £/ is the modal damping ratio. For certain 
forms of external loading, such as g(t) represented as 
a piecewise linear function of time, an exact explicit 
solution is available. This approach is described in 
reference 16 and is used in the NASTRAN® computer 
program (ref. 17). 

Equations (2.1) are the result of a given finite el- 
ement approximation designed to model the behav- 
ior of the dynamic system. The goal of the reduc- 
tion methods discussed in this section is to achieve 
an accurate approximation to the solution of equa- 
tions (2.1) with a small number of basis vectors. As 
discussed, the vibration modes are the most com- 
monly used basis functions in linear structural dy- 
namics. There are two cases, however, where a large 
number of modes are required for an accurate solu- 
tion of equations (2.1), and therefore the performance 
of the mode displacement method is poor. In the first 
case, if the structure is loaded in an impulsive man- 
ner, many high-frequency modes tend to be excited. 
These high-frequency modes must be included in the 
analysis since their contribution to the total response 
is significant. In the second case, if the response of 
the structure contains a large static component, the 
linear combination of vibration modes can do a poor 
job of approximating the static deflection shape. The 
reduction methods discussed in sections 2.2.2, 2.2.3, 
and 2.2.4 alleviate this second accuracy problem with 
the mode displacement method. 

2.2.2. Mode Acceleration Method 

To alleviate the poor accuracy of the mode dis- 
placement method due to its poor representation 
of the static component in the response, a method 
was proposed by Williams and Jones (ref. 18) called 
the mode acceleration method, which is described 
in its modern computational forms in references 16 


and 19. The mode acceleration method can be de- 
rived by rewriting equations (2.1) as 

u(t) = $(i)K -1 f-K -1 Cu-K _1 Mu (2.11) 

The first term in equations (2.11) is the quasi-static 
solution with load amplitude determined by g(t ). 
This term is calculated by solving the equations 

Ku s = f (2.12) 

This solution is carried out in the standard way by 
first factoring K into a product of upper and lower 
triangular matrices and then performing a forward 
and backward substitution operation to obtain u s . 
The other two terms are calculated with the solution 
for u and ii from the mode displacement solution. 
In these terms K -1 is calculated as follows. Equa- 
tions (2.3) can be rewritten in matrix form as 

2 = M<£ (2.13) 

where here $ is the full set of rig eigenvectors and 

fi -2 is 




n 


-2 


U/r) 



(2.14) 


With the eigenvectors scaled so that 

$ r M<F - I (2.15) 

equation (2.13) can be written as 

$ T K$rr 2 = i (2.16) 

Premultiplying by (<F r ) _1 and postmultiplying by 
yields 

K$ft- 2 $ T = I (2.17) 

or 

K -1 = (2.18) 

When contains less than the full n g eigenvectors, 
this expression for K -1 is only approximate. How- 

ever, since ii is obtained from the mode displacement 
solution based on n r modes (<&q), K -1 Mii is exactly 
equal to and no approximation results from 

introducing equations (2.18). For the damping term, 
introducing equations (2.18) with n r vectors in <I> is 
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not exact. However, this is a convenient approxima- 
tion especially when modal damping is used. Consis- 
tent with these considerations, equations (2.11) can 
be rewritten as 

u(t) = 0(t)K -1 f - <f>n~ 2 Cq - $rr 2 q (2.19) 

The key to the effectiveness of this method is that 
the static solution is included explicitly in the solu- 
tion. It is also simple to apply, since it essentially 
just superimposes the static and mode displacement 
solutions. Since q and q are obtained from the mode 
displacement solution, u and u are identical to the 
values obtained in the mode displacement method. 

2.2.3. Static Mode Method 

An alternative approach to the mode acceleration 
method that accounts for the static solution slightly 
differently is termed the static mode method herein. 
In this method, the static solution is included as 
an additional “mode” in forming the reduced equa- 
tions (2.5). The procedure begins by calculating a set 
of n r — 1 eigenvectors with equations (2.3). Then 
the static solution is calculated as 

K3>i = f ( 2 . 20 ) 

To improve the orthogonality of the basis vectors, 
the components of the vibration modes are removed 
from the static solution by use of the Gram-Schmidt 

process _ _ 

(j) x = ^ - 4>c (2.21) 

where ^ 

c = <£ T M0! (2.22) 

The vector ~4> x is then concatenated with to yield a 

new $ which is the complete basis. Equations (2.5) 
now become coupled and can be solved directly or 
reduced to an uncoupled form with the following 
procedure. First the reduced eigenvalue problem 

MZA + KZ-0 (2.23) 

is solved for the n r x n r diagonal matrix of eigen- 
values A and the n r x n r matrix of eigenvectors Z. 
Now a new set of basis vectors can be written as 

$ = $Z (2.24) 

When is substituted for <I> in equations (2.6), (2.7), 

(2.8), and (2.9), an uncoupled system results when C 
is of the special form described in section 2.2.1. 

The static mode method is similar to the mode 
acceleration method in that the static displacement 


vector is explicitly included in the solution. However, 
in the mode acceleration method the amplitude of the 
static displacement vector is not an unknown but 
is determined by g(t ), whereas in the static mode 
method the amplitude varies to possibly improve 
the solution. Also, this static displacement vector 
participates in the calculation of u and ii to possibly 
improve them as well. 

2.2.4. Ritz- Wilson- Lanczos Method 

A fourth method, which has become popular in 
the past few years is termed the Ritz-Wilson-Lanczos 
(RWL) method and is described in references 20, 21, 
and 22. Instead of using eigenvectors of the structure, 
this method uses a set of Lanczos vectors to form the 
reduced equations. The algorithm used here follows 
that in reference 20. The first vector is obtained by 
solving the static equations (2.20) and then scaling 
so that ^ 



The vectors i = 2,...,m arc obtained as follows. 
First 

K4>i = M<t>i-\ (2-26) 

is solved for (f )( . Then <pj is made NTorthogonal with 
respect to all previously generated vectors by using 
a Gram-Schmidt process 


i—l 


Gj 4>j 

(2.27) 

j= 1 


where ^ 


Gj — $ j 

(2.28) 

and scaling gives 


^ >t (~T ~ \l/2 

(2.29) 

(4>i M0i) 



It has been pointed out in many references (e.g., 
ref. 21) that the M-orthogonalization (eq. (2.27)) is 
theoretically required only with respect to the two 
previously computed vectors. However, it is also 
well-known that round-off errors cause the Lanczos 
vectors to become less and less orthogonal. Per- 
forming the Gram-Schmidt operation with respect to 
all previously generated vectors will not ensure the 
M-orthogonality of the vectors. However, it can im- 
prove the orthogonality in some cases. 

Following reference 20, a final step is performed in 
generating the basis vectors to produce an uncoupled 
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dynamic system. Of course this_is useful only when 
the system is undamped or C is assumed to be 
diagonal. As in the static mode method, a reduced- 
order eigenvalue problem is solved (eqs. (2.23)), and 
a new set of basis vectors produced. The process of 
explicitly computing the reduced stiffness and mass 
matrices required in equations (2.23) helps alleviate 
the problems caused by the lack of orthogonality of 
the Lanczos vectors. The matrices M and K in 
equations (2.23) are assumed to be full. That is, 
no assumptions arc made that particular terms in M 
and K are zero based on the properties of the vectors. 

2.3. Transient Response Solution Method 


When the reduction methods are used and general 
damping is included in the model, equations (2.5) are 
coupled. In principle any of the implicit or explicit 
numerical integration methods used for solving equa- 
tions (2.1) could be used to solve equations (2.5). In 
contrast to equations (2.1), however, equations (2.5) 
are low order, not stiff, and the primary concern is 
accurately integrating every equation in the system. 
Therefore an integration method which reduces trun- 
cation errors in the solution is highly desirable. Ac- 
curacy is especially important in sensitivity analy- 
ses because errors in the solution process are usually 
magnified in the calculation of derivatives. 

An approach that allows the use of moderately 
large time steps and makes the truncation error very 
small is called the matrix series expansion method in 
reference 23 and the transfer matrix method in refer- 
ences 24 and 25 when applied to structural dynamics 
problems and it is often referred to as a Taylor se- 
ries method in numerical analysis texts (e.g., ref. 26). 
This method expands the solution in a Taylor series 
where the number of terms determines the accuracy 
of the approximation. With this series, an expression 
can be written for the solution at time t + At in terms 
of the solution and load at time t as follows: 


fq(t + At)l fW n W 12 l/q(t)\ 

\q(f + At)/ [W 21 W 22 J\q(t)J 


Nn Nj 2 
N 2 i N 22 


0(*)f\ 

m f/ 


(2.30) 


It has been assumed here that the time variation of 
the load g(t) is approximated as a piecewise linear 
function of time, and therefore the second and higher 
order derivatives equal zero. Expressions for W x j 
and N,j are fairly complex and can be found in 
references 23 and 24. The values of the coefficients 
W ij and N,j depend on the number of terms taken 
in the series. 


The convergence properties of the W ^ series for 
an undamped, single-degree-of-freedom case can be 
studied by considering the following Taylor series 
expansion: 


A , . (wA t) 2 (ujAt) 4 

cos ijjAt = 1 - — + ^ 

2 4! 

, (wAt£ 

8 ! 


{uAtf 

6! 


(2.31) 


It is well-known that round-off errors due to finite 
precision arithmetic will cause large errors in this 
series for “large” values of coAt. Thus if a; is taken as 
w 7 j r , an upper bound on At can be estimated based 
on round-off error. In practice At usually needs 
to be much smaller than this upper bound value 
for two reasons. The first reason is that the input 
load history may be a complicated function of time, 
and At must be small enough to accurately sample 
this loading. The second, more important, reason is 
that At must be small enough to accurately sample 
the history of the output quantities. If At is larger 
than the smallest significant period of response, peak 
values of the response quantities will likely be missed. 
Accordingly, in the studies reported herein At was 
taken to be approximately one-eighth of the smallest 
period. Since the number of terms in the series has 
only a very small effect on the computational cost 
of the method, 50 terms were used in this study to 
make the truncations errors negligible. 
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Chapter 3 

Critical Point Constraint 


3.1. Constraint Formulation 

The general form of the constraint equation is 

ffi(x,t)<0 (3.1) 

where x is a vector of design variables and t is time. 
An effective approach for ensuring that this con- 
straint is satisfied for all values of t is the “critical 
point constraint” approach described in reference 27, 
pages 168-169. In this approach a set of peak val- 
ues of the function gi (denoted critical points) is se- 
lected. An obvious point to include is the time with 
the “worst” value of g t . However, if only this point is 
included, an optimization process modifying a struc- 
ture based on this information might unknowingly 
produce a design where the constraint is violated 
at another time point. To guard against this pos- 
sibility, a number of important peaks are selected. 
References 28 and 29 consider in detail the efficient 
location of critical points in large-scale structures 
problems with many constraints. This chapter 
presents a method for selecting the most important 
peaks as critical points. 

In the work reported herein, constraints are as- 
sumed to be placed on the displacements, velocities, 
accelerations, and stresses in the structure. All these 
constraints are treated similarly. Thus the critical 
point constraint formulation is illustrated for the case 
of displacements. Constraints are placed on selected 
displacements such that 

Mx,£)| < ^ allow ( 3 - 2 ) 

where u z are the displacements at specific points in 
the structure and u a u ow IS th e absolute maximum 
allowable value of the displacement. The critical 
values of this constraint occur at points in time where 
U{ has the largest magnitude. These are identified 
by examining every value of u along the response 
history. In the implementation here, each constraint 
is assumed to have a specified number of critical 
points; five critical points for each U{ are selected. 
Values of u where dujdt = 0 or values of u at the 


end points of the time interval are local maxima of 
g i and are termed candidate critical points. 

3.2. Selection of Critical Points 

The procedure for selecting the critical points 
from these candidates can best be explained by refer- 
ring to an example displacement time history shown 
in figure 3.1. The critical points are labeled with 
numbers and a few of the many candidate critical 
points are labeled with letters. The selection crite- 
ria applied to every candidate critical point are ex- 
plained as follows by considering these few candidate 
points. Candidate critical points a and c were dis- 
carded because the absolute values of the displace- 
ments at these points were smaller than those at 
the five other critical points. The criterion for dis- 
carding candidate points b, d, and e is slightly more 
complicated. From figure 3.1 it can be seen that all 
three candidate points have larger displacement mag- 
nitudes than that of critical point 1, for example. 
However, candidate points b, d, and e are all part 
of “major” peaks where a critical point is selected. 
A second criterion applied to the selection process 
is a requirement that only one critical point from 
each major peak be selected. This ensures that the 
critical points represent the total dynamic response 
rather than just the high-frequency undulations on, 
at worst, a single major peak. 



Figure 3.1. Example displacement time history illustrating 
critical point constraint selection process. 

A major peak is identified with the following 
procedure. Whenever a critical point is selected after 
comparing its magnitude with that at other critical 
points, a special screening process is activated. This 
screening process tests the displacement at every 
subsequent time point to determine if it differs from 
that at this last selected critical point by at least 
a specified percentage (25 percent for the studies 
reported herein). If so, all subsequent time points 
are no longer considered part of the current major 
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peak. Any candidate critical points identified while 
this special screening process is in effect are compared 
only against the last selected critical point. 

An example is the major peak in figure 3.1 which 
contains points d and 4. In the selection process, 
point d is initially selected as a critical point and 
the screening process is activated. The three points 
where du/dt = 0 between point d and point 4 are 
recognized to be part of the same major peak as d, 
but since the magnitude of the displacements at these 
points is smaller than at point d, they are discarded. 
Point 4 is also part of the same major peak as point d, 
but since the displacement magnitude there is larger 
than at point d, it replaces point d as a critical point. 
Before the next candidate critical point is considered, 
the displacement has changed from that at point 4 
by more than 25 percent and therefore is considered 
to be on a new major peak. 

3.3. Derivatives of Critical Point 

Constraints 

Once the critical points have been identified for 
the nominal design, these can be used in calculating 
sensitivities. Reference 27 demonstrates that the 
change in time location of critical points can be 
neglected in calculating derivatives of peak values 
with respect to design variables by examining the 
expression for the total derivative of g t with respect 
to a design variable x. Considering a constraint 
g(x,t) at a critical time t c gives 

dg(x,t c ) _dg dgdtc 

dx dx dt dx (6 ' 6) 

The last term in equation (3.3) is always zero be- 
cause at interior critical points dg/dt = 0 and on the 
boundary dt c /dx — 0. Accordingly, the sensitivity 
calculations need to be performed only at the spe- 
cific times where the critical points have been iden- 
tified. This can result in a considerable savings in 
computational time, especially when there are many 
constraints, many time points, or many basis func- 
tions used to represent the response. The details of 
each sensitivity calculation method are discussed in 
the next chapter. 
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Chapter 4 

Methods For Calculating Sensitivities 

4.1. Finite Difference Methods 

4.1.1. Forward and Central Difference Operators 

Both the forward difference and central difference 
methods have been used in this study to calculate 
sensitivities. The well-known forward difference ap- 
proximation to du/dx , 

A u u(x + Ax) — u(x) . 

Ax Ax 

and the central difference approximation, 

A u u(x + Ax) — u(x — Ax) 

Ax = 2Ax * j 

are used. The truncation error for the forward 
difference approximation is 

e r (Ax)= ^^(x + CAz) (0 < C < 1) (4-3) 
and for the central difference approximation is 

e r (A*)= (M!^( x + C Ax) (0 < < < 1) 

(4.4) 

In applying equations (4.1) and (4.2), the selec- 
tion of difference step size Ax is a concern. Selection 
of a large step size results in errors in the deriva- 
tive due to truncation of the operator (eqs. (4.3) and 
(4.4)). Selection of a small step size can lead to er- 
rors in the derivative due to the limited floating point 
precision of the computer or algorithmic inaccuracies 
in calculating u (condition errors). It is not uncom- 
mon with the forward difference method (eq. (4.1)) 
that no acceptable value exists for Ax to produce an 
accurate value of du/dx considering the conflicting 
requirements of minimizing truncation and condition 
errors. Because the truncation error associated with 
the operator of equation (4.2) is typically less than 
that of equation (4.1), it is possible to use a larger 


finite difference step size. The larger Ax reduces the 
condition error from the function evaluations and re- 
sults in a more accurate value of du/dx . However, 
the necessity of two function evaluations needed for 
equation (4.2) makes the procedure computationally 
more costly. 

4.1.2. Using Vibration Modes as Basis 

Functions 

For many of the studies herein the natural vibra- 
tion mode shapes are used as basis functions to rep- 
resent the transient response. In calculating the re- 
sponse of the perturbed design in equation (4.1) and 
the two perturbed designs in equation (4.2), some 
computational savings are possible relative to the 
computations for the initial design. 

If the mode shapes for the initial design are used 
to represent the perturbed design, the cost of resolv- 
ing the eigenvalue problem is eliminated. However, 
the reduced set of equations_foi^the perturbed system 
must still be formed and M, C, and K are now full. 
This coupled system is then solved with the matrix 
series expansion method described in section 2.3. 

If the updated mode shapes for t he perturbed de- 
sign are used in the analysis, many eigensolution pro- 
cedures, such as the subspace iteration used here, can 
begin with the mode shapes from the initial design 
as approximations. Since the perturbation in the de- 
sign is small, the subspace iteration procedure con- 
verges rapidly. However, at least one factorization of 
K is required. For large finite element models, this 
can be the largest part of the computational cost. 
For most of the studies in this paper, the forward 
difference method used the initial mode shapes to 
represent the perturbed design. Because the central 
difference method was used for reference values of 
derivatives, updated mode shapes were calculated for 
the two required perturbed designs. In both cases, 
because of the critical point constraint formulation, 
the transformation from modal coordinates to phys- 
ical coordinates (e.g., displacements, stresses) is per- 
formed only at the critical times instead of at all time 
points. 

4.2. Semianalytical Methods 

The direct method for sensitivity calculation is 
derived by differentiating equations (2.1). The 
derivation presented here follows that in reference 27, 
pages 169-171. After differentiating equations (2.1) 
with respect to the design variable x the result is 
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„du _du T du df , . dM .. 

M — + C— + K— =—g(t) - — u 
ax ax ax ax dx 


dC . __ dK 
dx dx 


(4.5) 


This system of differential equations of order n g could 
be solved directly for the sensitivities du/dx , dix/dx , 
and du/dx. However, just as for the response equa- 
tions, it is more efficient to consider a reduced form 
of the sensitivity equations which can be obtained 
by differentiating equations (2.5) with respect to x 
to yield 


d A + cp + Kp=f g (t)- d -P 

ax dx dx dx dx 


functions of member cross-sectional area. In these 
cases, there are no truncation errors and large finite 
difference step sizes can be used to reduce the condi- 
tion error and produce accurate derivatives. 

Calculation of the first term in equation (4.7) 
and the second and third terms in equation (4.8) 
depends on the particular choice of basis functions 4>. 
Considerable reduction in computational cost results 
if the vectors $ are taken to be independent of x that 
is fixed. Methods which are less costly than exact 
methods are also available to approximate d$>/dx. 
Two semianalytical procedures which address these 
concerns are discussed in sections 4.2.1 and 4.2.2. 

4-2.1. Fixed- Mode Semianalytical Formulation 


dC . dK , x 
^ q -^ q <46) 


The first step in forming this equation is the calcula- 
tion of the derivatives of f, M, C, and K (eqs. (2.9), 
(2.6), (2.7), and (2.8)) with respect to x. Using equa- 
tions (2.9) gives the derivative of f with respect to x 
as follows: 


df d& T c rdf 

= f + <P 

dx dx dx 


(4.7) 


The force f is frequently not a function of the design 
variables; this simplifies equation (4.7). Also, with 
equation (2.6), the derivative of M with respect to x 
can be written as 


dM 

dx 





d* T 

dx 


M$ + <£ 7 M 


d& 

dx 


(4.8) 


Similar expressions can be written for the derivatives 
of C and K. 

The derivative dM/dx in equation (4.8) (simi- 
larly for dC/dx and dK/dx) is, in general, difficult 
to calculate because the finite element model may 
be composed of diverse element types whose prop- 
erties are complicated functions of the design vari- 
ables x. For this reason, these derivatives are often 
replaced with finite difference approximations. This 
combination of analytical differentiation of the re- 
sponse equations with finite difference matrix deriva- 
tives is known as a semianalytical approach. The 
semianalytical methods presented herein for calculat- 
ing transient response quantities all use the forward 
difference operator to approximate dM/dx, dC/dx, 
dK/dx, and df/dx. For several important classes of 
design variables, however, M, C, and K are linear 
functions of x. For example, M and K in a finite 
element model composed of truss members are linear 


If the basis vectors are assumed not to be func- 
tions of the design variables x, d4>/dx equals 0. This 
significantly simplifies equations (4.7) and (4.8). Af- 
ter forming the derivatives of f, M, C, and K, the 
right-hand side of equations (4.6) can be formed us- 
ing Q? and q from the solution of equations (2.5). 
The matrix series expansion method ensures that ac- 
curate values of q, q, and q are available for this 
step. Equations (4.6) can then be integrated to yield 
dq/dx, dq/dx, and dq/dx. Herein this fixed-mode, 
semianalytical implementation of the direct method 
is called the semianalytical method. 

4-2.2. Variable-Mode Semianalytical 

Formulation 

If the basis functions are assumed to be functions 
of x, the calculation of d4>/dx either exactly or ap- 
proximately is required to form equations (4.7) and 
(4.8). Vibration modes are the most popular basis 
functions, and the calculation of their derivatives has 
been studied extensively. Reference 30, for exam- 
ple, surveys several methods for calculating deriva- 
tives of vibration mode shapes from a computational 
point of view. One of the most popular methods, 
Nelson’s method (ref. 31), requires a factorization of 
the system equations for each eigenvector considered. 
This can be a considerable computational burden for 
large systems. Since the overall objective here is 
the accurate calculation of transient response sen- 
sitivities, not eigenvector sensitivities, it seems desir- 
able to investigate cheaper approximate methods for 
calculating d4>/dx. 

One approximate method for calculating the 
eigenvector derivatives is similar to the modal ap- 
proach for transient analysis. This modal method 
approximates each eigenvector derivative as a 
linear combination of the modes themselves. 
In many cases, however, a very large number of 
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eigenvectors are required for accurate derivatives. 
Furthermore, the eigenvector derivative approxima- 
tion produced by this method cannot improve the 
transient response sensitivities because they are 
contained entirely within the span of the modes 
themselves. 

A method proposed by Wang (ref. 32) to alle- 
viate the poor performance of the modal method 
also improves the transient sensitivities. This modi- 
fied modal method is derived by first differentiat ing 
equations (2.3) to yield 


(K - uj M) 


,2i 

dx 


d w? 

— - 
dx 


dK 

dx 


+ Uj 


d M 


dx J 
(4.9) 

This equation cannot be solved directly since the left- 
hand side is singular. Wang’s approach, however, 
was to calculate a pseudostatic solution to this equa- 
tion by neglecting M on the left-hand side of equa- 
tions (4.9). The solution to this pseudostatic equa- 
tion introduces the change in basis associated with 
changes in the design variables and is significant in 
improving the transient response sensitivities. The 
mode shape derivative can then be written as 


with similar expressions for du/dx and du/dx. The 
calculation of stresses begins with 

<j = Su (4.14) 


where S is the stress transformation matrix. Substi- 
tuting equation (2.4) yields 


(j — S<t>q 


(4.15) 


when equation (4.15) is differentiated, the stress 
sensitivities can be written as 


dcr 

dx 


_ ^ dq dS 

+ —$q + 
dx dx 



(4.16) 


The matrix dS/dx is approximated with the forward 
difference operator. Because of the critical point con- 
straint formulation, the transformation from these 
modal quantities to physical displacements, veloci- 
ties, accelerations, and stresses is performed only at 
the critical times. 


4.2.4 ■ Mode Acceleration Method 


d& 3 

dx 


d*j 

dx 


+ 51 A i k 

1 




(4.10) 


where ( d&j/dx) s is the pseudostatic contribution. 
The coefficients A jk are obtained by substituting 
equation (4.10) into equations (4.9), multiplying by 
4>J, and simplifying as 


The mode acceleration method was presented in 
chapter 2 as a technique for improving the dynamic 
displacements and stresses when the static compo- 
nent is significant. It is also possible that it can im- 
prove the sensitivities of displacements and stresses. 
An expression for the sensitivities using the mode 
acceleration method is obtained by first rearranging 
equations (4.5) to yield 


A i k — 


4*1 { 


dK 

dx 


, .2 dM 
“j HP 


) 




(Mi) (4.H) 


or 


du(t) 

dx 


K 


-l 


df . . dK dii dC . 

dx 9 ^ dx U dx dx U 


, , dii dM .. 

*r u 


(4.17) 


= ~\*Jl£*i (* = 7 ) < 4 ' 12 > 

Given these approximate values of eigenvector deriva- 
tives, equations (4.7) and (4.8) can be formed. Then 
equations (4.6) can be solved for dq/dx, dq/dx, 
and dq/dx just as in the fixed-mode, semianalytical 
method. 


4-2.3. Recovery of Physical Sensitivities 


Given dq/dx, the derivative of the physical dis- 
placement vector du/dx can be written as 


du 

dx 




(4.13) 


If a reduced basis approximation is applied uniformly 
to every term in equations (4.17), the resulting du/dx 
would agree with that obtained from the solution 
of equations (4.6). The objective of a mode accel- 
eration solution is to selectively apply the modal 
approximation to equations (4.17) with the goal of 
improving the values of du/dx. In applying the 
mode acceleration method to the transient response 
problem (eqs. (2.11)), u and u are obtained from 
the mode displacement method. Here, in apply- 
ing the mode acceleration method to the sensitivity 
equations, u is obtained from the solution to equa- 
tions (2.19) and du/dx and dii/dx are obtained from 
the solution to the mode displacement, semianalyti- 
cal equation (eqs. (4.6)). In the derivation here, the 
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modes are assumed to be fixed. Substituting equa- 
tions (2.19) and u = 4>q, ii = 3>q, du/dx = $dq/dx, 
and du/dx = &dq/dx into equations (4.17) yields 


du(f) 

dx 


= K ] 


di dK 

Tx s(t) ~ ta 


'g(t) K~U 


-W 2 Cq- W z q -C$ 


2” 


( dq 

dx 


dc _. A/rA dq dM 

4>q - Mf>- — 3>q 

dx dx dx 


(4.18) 


method is that it is not as effective as the method 
of equations (4.19). The mode acceleration method 
in the overall finite difference procedure provides a 
good approximation to the first term (pseudostatic 
term) in equations (4.19). However, the key effects 
of the K^ 1 (dK/dx)^n~ 2 and the K- l (dM/dx)<$> 
terms are completely neglected. 


The modal approximation for K _1 (eqs. (2.18)) is 
introduced into all terms in equation (4.18) that 
involve damping just as in the mode acceleration 
solution described in chapter 2. It was also pointed 
out in chapter 2 that K _1 M<£ in equation (4.18) 
is exactly Based on these considerations, 

equation (4.18) can be simplified to yield the mode 
acceleration solution of the sensitivity equations 


t/u(£) 

dx 


K 


df _4 K k -i f 

dx dx 


9(t) 




-2 


dx dx 


+ K 


dK 

dx 




dM 
* 

dx 


q-$fr 2 c^ 

dx 

q 

dx 


(4.19) 


The key to the effectiveness of this mode acceleration 
sensitivity method is the usage of the exact K -1 
in the calculation of the K -1 (dK/dx)&Q~~ 2 and 
K _1 (dM/da:)* terms. The explicit calculation of 
these terms expands the basis beyond the span of 
the modes in a manner similar to the pseudostatic 
term in the modified modal method described in the 
previous section. When equations (4.19) are used, 
the stress sensitivities can be calculated as 


da gdu dS 

dx dx dx 


(4.20) 


where u is obtained from equations (2.11). 

It is worthwhile to contrast the sensitivity ap- 
proach of equations (4.19) with an alternate approach 
of a fixed-mode, overall forward difference method 
with the response quantities calculated with a mode 
acceleration method. This overall forward difference 
approach has one obvious drawback. The mode ac- 
celeration method requires the costly factorization of 
K for the perturbed design. Therefore, much of the 
cost savings achieved by keeping the modes fixed is 
lost. A second defect of the overall forward difference 
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Chapter 5 

Numerical Studies 

The different transient response methods de- 
scribed in chapter 2 and the sensitivity calculation 
methods described in chapter 4 are applied to three 
example problems in this chapter. The three example 
problems are small but they all have certain charac- 
teristics which complicate the dynamics and sensi- 
tivity calculations. The first example is the five-span 
beam with relatively closely spaced frequencies and 
loaded with a moment applied at a single point. As 
a result, many modes participate in the dynamic re- 
sponse. The second example is the delta wing loaded 
with a uniform pressure load. Although the higher 
frequency modes are not significantly excited by this 
loading, the analysis is complicated by the laminated 
plate elements in the model and the sensitivity anal- 
ysis is complicated by the lamina thickness design 
variables considered. The third example is a can- 
tilever beam with a stepped cross section loaded with 
an applied rotation at the root. This loading is iner- 
tial, depends on the mass, and therefore also depends 
on the values of the design variables. The first two 
examples consider point mass and standard thick- 
ness design variables. The cantilever beam example 
also includes so-called shape design variables (sec- 
tion lengths) that are known to cause difficulties in 
the sensitivity analysis in some cases. 

One of the key questions addressed in this chap- 
ter is how well a particular set of basis vectors rep- 
resents the full system of order n g . This full system, 
however, is the result of a particular finite element 
discretization. Thus the accuracy of the response or 
sensitivities as a function of the finite element mesh is 
also an appropriate question. This question is espe- 
cially important when a large number of basis vectors 
(n r close to n g ) are required for an accurate solution 
in a problem with a given finite element mesh. Either 
the basis vectors are doing a very poor job of span- 
ning the solution space or the loading is legitimately 
exciting this high-frequency behavior. 

In this chapter two terms are used to describe 
these studies which consider the dynamic response 
as a function of the number of basis vectors or the 
number of finite elements in the model. The effect 
of the number of basis vectors on the accuracy of 


the response or sensitivities for a given finite element 
mesh is called a “modal convergence” study. In this 
case, the goal is for the n r basis vectors to provide 
an accurate solution to the approximate equations of 
order n g . The question of whether the finite element 
model associated with this system is an accurate 
representation of the continuum is addressed in a 
“mesh convergence” study. In some cases, it will be 
shown that the modal convergence is strongly related 
to the mesh convergence. That is, when a large 
number of basis vectors are required for an accurate 
solution for a given finite element mesh, the finite 
element mesh is doing a poor job of representing the 
continuum. In other cases, even though the finite 
element mesh is providing an accurate representation 
of the continuum, some sets of basis vectors are doing 
a poor job of representing the response or sensitivities 
for this n 9 th order system. 

Several additional comments on the concept of 
a modal convergence study are in order. Clearly 
the use of the term convergence is imprecise because 
the accuracy of the approximate solution with dif- 
ferent numbers of modes is compared only with the 
fi n ite-degree-of- freedom solution rather than the con- 
tinuum solution. However, it is assumed that an 
“acceptable” finite element model must do a good 
job of representing the low-frequency modes of the 
structure. Therefore, the accuracy of the dynamic 
solution with a small number of modes from the 
finite-degree-of-freedom model is a very reasonable 
approximation to the accuracy of the dynamic solu- 
tion with a small number of modes calculated from a 
continuum model. Thus the convergence of the solu- 
tion as a function of the number of modes calculated 
from the finite element model is a reasonable approxi- 
mation to the true modal convergence obtained when 
the modes are calculated from the continuum model 
as long as the number of modes considered is small. 
Furthermore, if the number of modes required for an 
accurate calculation of either the response or sensi- 
tivities is not small, the basis vectors or the method 
will be considered poor. 

5.1. Five-Span Beam Example 

The first example considered is a five-span, planar 
beam example taken from reference 33 and shown 
in figure 5.1. An initial investigation of the dis- 
placement transient response of this problem was 
also considered in reference 34. In most of the 
studies, the beam is modeled with three beam fi- 
nite elements per span resulting in 26 unconstrained 
degrees of freedom. The effect of finite element 
discretization is considered by developing alternate 
models with 6, 9, 12, ... elements per span. As 
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shown in figure 5.1, translational and rotational vis- 
cous dampers were also added to the beam. These 
devices are representative of velocity feedback con- 
trollers which might be added to flexible structures. 
Cases with and without dampers were considered. 
The numerical values of the damping coefficients 
from reference 33 of c\ = 0.008 sec-lb/in. and c<± — 
1.2 sec-lb were used. In one example, modal damp- 
ing with = 0.005, which is intended to represent 
typical structural damping, was used instead of the 
discrete dampers. A case was also considered where a 
1.0-lb mass (approximately 20 percent of the mass of 
the beam) was added to the beam at the location of 
the translational damper. The eigenvalues for three 
cases using the model with three elements per span 
are shown in table 5.1. The additional point mass 
has a significant effect on the frequencies, whereas 
the dampers have little effect. The effect on frequen- 
cies of increasing the number of elements per span in 
the finite element model is shown in table 5.2. It can 
be seen that the lowest 10 frequencies are fairly well 
converged even for the model with three elements 
per span. In the transient analysis, the applied load- 
ing for all problems consisted of a point moment of 
0.04405 in-lb applied at the right end of the beam. 
Two different time functions for this load, a step and 
a ramp (shown in fig. 5.1), were considered. 
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Figure 5.1. Five-span beam with applied end moment. 


5.1.1. Beam Dynamic Response 

The first part of this study focused on the tran- 
sient response of the beam with the mode dis- 
placement, mode acceleration, static mode, and 
Ritz-Wilson-Lanczos (RWL) methods. Displace- 
ment, velocity, acceleration, and stress resultant re- 
sponse quantities are considered. For this beam ex- 
ample, all these response quantities are taken at a 
location 10.0 inches from the left end of each span. 
This point is the end of the first element in each span 
when three elements per span are used in the model. 

5. 1.1.1. Character of response. In the first 
case, the ramp loading was applied to the undamped 


Table 5.1. Eigenvalues For Three Five- Span Beam Cases 



Undamped 

Damped with 
point dampers 

Undamped 
with point 
mass 

Mode 

Frequency, Hz 

Frequency, Hz 

Damping ratio 

Frequency, Hz 

1 

1.1707 

1.2210 

0.0851 

0.9401 

2 

1.2991 

1.2926 

.0352 

1.2594 

3 

1.6254 

1.6298 

.0690 

1.5445 

4 

2.0491 

2.0910 

.0590 

1.8005 

5 

2.4628 

2.5497 

.0958 

2.3729 

6 

4.7343 

4.8426 

.0044 

4.2327 

7 

5.0105 

4.9785 

.0413 

4.8858 

8 

5.6472 

5.7703 

.0126 

5.6400 

9 

6.4153 

6.4178 

.0407 

5.9261 

10 

7.1274 

7.2229 

.0193 

6.8762 


Table 5.2. Beam Frequencies With Different Numbers of 
Finite Elements Per Span 



Frequency, Hz, for — 

Mode 

3 elements 

6 elements 

9 elements 

12 elements 

1 

1.1707 

1.1698 

1.1698 

1.1698 

2 

1.2991 

1.2979 

1.2978 

1.2978 

3 

1.6254 

1.6230 

1.6229 

1.6229 

4 

2.0491 

2.0445 

2.0442 

2.0442 

5 

2.4628 

2.4547 

2.4542 

2.4542 

6 

4.7343 

4.6828 

4.6798 

4.6792 

7 

5.0105 

4.9504 

4.9469 

4.9462 

8 

5.6472 

5.5652 

5.5601 

5.5593 

9 

6.4153 

6.3053 

6.2980 

6.2967 

10 

7.1274 

6.9974 

6.9874 

6.9857 


beam modeled with three elements per span. All 
26 modes were used in the analysis. Time histo- 
ries of selected displacement, velocity, acceleration, 
and bending moment components are shown in fig- 
ures 5.2, 5.3, 5.4, and 5.5, respectively. The displace- 
ment history (fig. 5.2) is relatively smooth indicat- 
ing that only the low-frequency modes of the beam 
are contributing to the response. The velocity and 
bending moment response histories are more jagged 
indicating participation by higher frequency modes. 
The acceleration history (fig. 5.4) is extremely jagged 
with contributions from the highest frequency modes 
represented by the finite element model. 

The impulsive nature of the step load makes 
the higher frequency modes much more important. 
This can be seen in figure 5.6 where the time his- 
tory of velocity in the second span (fig. 5.1) U 2 is 
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shown. By comparing this velocity history with that 
in figure 5.3 the increased importance of the high- 
frequency modes becomes obvious. 



Time, sec 

Figure 5.2. Time history of displacement U 2 for five-span 
beam subjected to transient end moment. Ramp 
load; undamped beam. 



Figure 5.3. Time history of velocity ii 2 for five-span beam 
subjected to transient end moment. Ramp load; 
undamped beam. 

The addition of the point dampers shown in fig- 
ure 5.1, on the other hand, tends to reduce the im- 
portance of the high-frequency modes. This is shown 
in figure 5.7 where again ii2 is shown. Comparing fig- 
ures 5.7 and 5.3 shows that the velocity history for 
the damped case is significantly smoother than for 
the undamped case. 

This changing character of the time histories with 
temporal or spatial differentiation of the response 
function or the addition of dampers is expected. The 
implications of this phenomenon on calculating the 
sensitivities of these response quantities are discussed 
below. 

5. 1.1. 2. Modal convergence. When vibration 
modes or other functions are used to reduce the ba- 
sis in a transient response problem (eq. (2.4)), the 



Figure 5.4. Time history of acceleration 112 for five-span 


beam subjected to transient end moment. Ramp 
load; undamped beam. 



Figure 5.5. Time history of bending moment in span 5 for 


five-span beam subjected to transient end mo- 
ment. Ramp load; undamped beam. 



Figure 5.6. Time history of velocity U 2 for five-span beam 
subjected to transient end moment. Step load; 
undamped beam. 

key question is how many modes are required for 
an accurate solution. This section addresses that 
question for the five-span beam example with the 
response calculated with mode displacement, mode 
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acceleration, static mode, and RWL methods. Un- 
less otherwise stated, all the response quantities are 
considered at critical times selected with the methods 
discussed in chapter 3. 



Figure 5.7. Time history of velocity ii ‘2 for five-span beam 
subjected to transient end moment. Ramp load: 
damped beam. 
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Figure 5.8. Modal convergence of selected displacements for 
five-span beam. Ramp load; undamped beam; 
mode displacement method. 

The baseline case of ramp loading applied to the 
undamped beam modeled with three elements per 
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Figure 5.9. Modal convergence of selected velocities for five- 
span beam. Ramp load; undamped beam; mode 
displacement method. 

span is considered first. Figure 5.8 shows the conver- 
gence of selected displacements at critical points as 
a function of number of modes. The displacement 
critical point combinations were selected to be rep- 
resentative of both the largest and smallest critical 
values. In figure 5.8 and in all the other figures show- 
ing convergence of response quantities or sensitivi- 
ties, the figure key indicates the quantity and the 
time of occurrence in seconds. In all cases the con- 
vergence is very good with approximately 10 modes 
yielding a converged solution. Figure 5.9 shows a 
similar plot for velocities. Again the convergence is 
good. The modal convergence for accelerations, how- 
ever, is poor as shown in figure 5.10. Figures 5.11 and 
5.12 show the modal convergence of selected bending 
moments and shear forces, respectively, and again, 
the convergence is poor. 

To possibly alleviate this poor convergence, the 
alternate reduction methods discussed in chapter 2 
were applied to this problem. The modal conver- 
gence for displacements calculated with the mode 
acceleration method (fig. 5.13) is even better than 
that found with the mode displacement method. The 
convergence of bending moments and shear forces 
has improved dramatically from the mode displace- 
ment results as can be seen in figures 5.14 and 5.15. 
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Figure 5.10. Modal convergence of selected accelerations for 
five-span beam. Ramp load; undamped beam; 
mode displacement method. 
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Figure 5.11. Modal convergence of selected bending mo- 
ments for five-span beam. Ramp load; un- 
damped beam; mode displacement method. 
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Figure 5.12. Modal convergence of selected shear forces for 
five-span beam. Ramp load; undamped beam; 
mode displacement method. 
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Figure 5.13. Modal convergence of selected displacements 
for five-span beam. Ramp load; undamped 
beam; mode acceleration method. 
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As mentioned in chapter 2, the mode acceleration 
method does not apply to the calculation of veloci- 
ties and accelerations. 

A similar improvement was noted with the static 
mode method. As an example, consider the excellent 
convergence of shear forces shown in figure 5.16. 
However, the addition of the static solution provides 
no improvement in the convergence of acceleration as 
shown in figure 5.17. 

The RWL method is attractive because of the 
significantly reduced cost of calculating the vectors 
compared with solving the eigenproblem. In this 
five-span beam example, the modal convergence is 
also as good as the mode acceleration or static mode 
methods. The good convergence of the shear forces 
is shown in figure 5.18. Like the other reduction 
methods, however, the convergence of accelerations 
is poor (fig. 5.19). 

The modal convergence of the response quanti- 
ties for the step-loaded case is generally much poorer 
than for the ramp-loaded case. The convergence of 
the displacements is reasonably good. Convergence 
of velocities, accelerations, and stresses, however, is 
poor. This poor convergence is not surprising con- 
sidering the “jaggedness” of the velocity time his- 
tory shown in figure 5.6. As an example, two figures 
with the convergence of bending moments plotted as 
a function of the number of modes are shown. The 
first, figure 5.20, shows the bending moments cal- 
culated with the mode displacement method. Con- 
vergence is poor but this is not surprising since the 
convergence was poor with the mode displacement 
method for the ramp-loaded case (fig. 5.11). For 
this ramp-loaded case, the convergence of the bend- 
ing moments improved dramatically when the mode 
acceleration method was used as can be seen in 
figure 5.14. Although convergence is improved for 
the step-loaded case by using the mode acceleration 
method (fig. 5.21), the convergence is still fairly poor. 

Judging from the velocity time history in fig- 
ure 5.7, it might be expected that including damp- 
ing would improve the modal convergence of the re- 
sponse quantities. For the ramp-loaded, undamped 
case, the poorest convergence was for the acceler- 
ations (fig. 5.10). For the ramp-loaded case with 
discrete dampers, there is an improvement in modal 
convergence as seen in figure 5.22. For the case with 
0.5 percent modal damping there is also a slight im- 
provement in modal convergence. However, in nei- 
ther case does the damping completely alleviate the 
poor convergence. 

All the previous convergence results are at critical 
points located by the method described in chapter 3. 
When a different number of modes is used in the anal- 
ysis, the critical time for a particular critical point 
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Figure 5.14. Modal convergence of selected bending mo- 
ments for five-span beam. Ramp load; un- 
damped beam; mode acceleration method. 
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Figure 5.15. Modal convergence of selected shear forces for 
five-span beam. Ramp load; undamped beam; 
mode acceleration method. 

usually shifts slightly. Consequently, the results for 
a given response quantity-critical point combination 
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Figure 5.16. Modal convergence of selected shear forces for 
five-span beam. Ramp load; undamped beam; 
static mode method. 
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Figure 5.17. Modal convergence of selected accelerations for 
five-span beam. Ramp load; undamped beam; 
static mode method. 
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Figure 5.18. Modal convergence of selected shear forces for 
five-span beam. Ramp load; undamped beam; 
Ritz-Wilson-Lanczos method. 
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Figure 5.19. Modal convergence of selected accelerations for 
five-span beam. Ramp load; undamped beam; 
Ritz-Wilson-Lanczos method. 
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Figure 5.20. Modal convergence of selected bending mo- 
ments for five-span beam. Step load; un- 
damped beam; mode displacement method. 
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Figure 5.21. Modal convergence of selected bending mo- 
ments for five-span beam. Step load; un- 
damped beam; mode acceleration method. 


occur at different times depending on the number 
of modes used in the analysis (the values shown in 
parentheses in the figures are for the most refined so- 
lution). It is natural to ask whether response quanti- 
ties at fixed times plotted as a function of number of 
modes would show similar convergence. Figures 5.23 
and 5.24 show the modal convergence of selected ve- 
locities and bending moments, respectively, at fixed 
times. The mode displacement method was used in 
the analyses. The particular response quantities and 
times were selected to span the range between largest 
positive and negative values. As can be seen in fig- 
ure 5.23, the convergence of velocities is good. From 
figure 5.24, it can be seen that the convergence of the 
bending moments is poor and remarkably similar to 
the critical point convergence results (fig. 5.11). Thus 
it would appear that the critical point constraint for- 
mulation does not significantly affect the modal con- 
vergence of the response. 

5. 1.1.3. Mesh convergence. Table 5.2 shows 
the convergence of the lowest 10 frequencies as a 
function of the number of elements used to model 
each span of the beam. The convergence of these 
lower frequencies is rapid. The convergence of various 
response quantities as a function of the mesh is also 
a concern. 

The modal convergence cannot be uncoupled 
from the mesh convergence. This was discussed in 
reference 33 for the derivatives of damping ratios 
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Figure 5.22. Modal convergence of selected accelerations for 
five-span beam. Ramp load; discretely damped 
beam; mode displacement method. 
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Figure 5.23. Modal convergence of selected velocities for 
five-span beam at fixed time points. Ramp 
load; undamped beam; mode displacement 
method. 
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Figure 5.24. Modal convergence of selected bending mo- 
ments for five-span beam at fixed time points. 
Ramp load; undamped beam; mode displace- 
ment method. 


calculated with undamped vibration modes. For sev- 
eral cases, the modal convergence of the derivatives 
was poor. As the mesh was refined, convergence was 
achieved only when almost all the available modes 
were used in calculating the damping ratio. Clearly, 
this is an example where the modal basis provides a 
very poor approximation to the actual solution. 
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Figure 5.25. Modal convergence of selected shear forces for 
five-span beam modeled with six elements per 
span. Ramp load; undamped beam; mode 
displacement method. 

Figure 5.25 shows the modal convergence of shear 
force for the five-span beam modeled with six ele- 
ments per span and the transient analysis performed 
with the mode displacement method. The conver- 
gence for this case is just as poor as for the case with 
three elements per span shown in figure 5.12. How- 
ever, a plot of convergence of this shear force as a 
function of the number of elements per span, when all 
modes are used in each analysis (figure 5.26), shows 
good convergence. Clearly, the convergence of shear 
forces for the ramp-loaded five-span beam is similar 
to that reported for derivatives of damping ratios in 
reference 33; the vibration modes are simply doing a 
poor job of representing the solution. 

Figure 5.27, which shows the convergence of ac- 
celerations for the step-loaded beam as a function 
of elements per span, indicates a different behavior, 
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Figure 5.26. Convergence of selected shear forces as function 
of number of finite elements per span for five- 
span beam. Ramp load; undamped beam; 
mode displacement method. 
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Figure 5.27. Convergence of selected accelerations as func- 
tion of number of finite elements per span for 
five-span beam. Step load; undamped beam; 
mode displacement method. 


however. Here, the very poor convergence is due 
to the higher frequency modes being excited by the 
step loading. As the mesh is refined, the number of 
high-frequency modes increases, and these continue 
to have a significant contribution to the acceleration. 

In evaluating the accuracy of the sensitivity cal- 
culation procedures in section 5.1.2, particular atten- 
tion must be paid to the convergence characteristics. 
Some convergence problems such as those caused by 
the use of vibration mode shapes can be improved by 
the use of alternate basis functions. However, other 
convergence problems, such as for the accelerations 
in the step-loaded case, are inherent in the problem 
definition. 

5.1.2 . Sensitivities of Beam Dynamic Response 

In the previous section, the transient response of 
the five-span beam was considered in detail. In this 
section, the calculation of sensitivities of displace- 
ments, velocities, accelerations, and stresses with re- 
spect to various design variables is considered. 

5. 1.2.1. Design variables. Two different classes 
of design variables were considered. The first design 
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Figure 5.28. Effect of finite difference step size on accuracy 
of displacement derivative approximation with 
respect to point mass design variable. Ramp 
load; undamped beam. 
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Finite difference step size, in. 

Figure 5.29. Effect of finite difference step size on accuracy 
of displacement derivative approximation with 
respect to thickness design variable. Ramp 
load; undamped beam. 

variable is a concentrated mass m (initially zero) at 
the location of the translational damper. This de- 
sign variable was also considered in reference 33. The 
derivatives of the system mass and stiffness matrices 
with respect to this design variable are constant and 
zero, respectively. As a consequence, the derivatives 
of the system matrices required in the semianalyti- 
cal methods can be calculated exactly by a simple 
forward difference operator. The beam thicknesses 
in each of the five spans were also design variables. 
Derivatives with respect to the five thickness design 
variables showed similar characteristics. Herein, re- 
sults for derivatives with respect to the thickness in 
the rightmost span along with derivatives with 
respect to the point mass m are presented. 

5. 1.2.2. Effect of finite difference step size. The 
methods described in chapter 4 for calculating sensi- 
tivities all rely on finite difference operators at some 
stage in the algorithm. The forward and central 
difference methods rely on the operators in equa- 
tions (4.1) and (4.2) to calculate derivatives of re- 
sponse quantities. In the semianalytical methods, the 
derivatives of the system matrices are calculated with 
the forward difference operator in equation (4.1). In 
all cases the finite difference step size must be se- 
lected so that the operator provides a reasonable ap- 
proximation to the derivative. If the step size is too 
large, the error due to truncating the series approx- 
imation of the derivative is large. If the step size 
is too small, the numerical condition error in per- 


forming the function evaluations (dynamic analyses) 
becomes large. 

To assess the effect of step size on the calcula- 
tion of sensitivities for the five-span beam, deriva- 
tives were calculated with the three methods with 
various step sizes. In this study the beam was un- 
damped and the ramp loading was applied. All 26 
vibration modes were included in the analysis. Fig- 
ure 5.28 shows the estimated derivative of displace- 
ment U 2 at critical point 5 with respect to the point 
mass design variable m as a function of step size. 
As mentioned, the derivatives of the system matri- 
ces with respect to this design variable can be calcu- 
lated exactly in the semianalytical method. As a re- 
sult, the derivative estimated with the semianalytical 
method is approximately constant for the six-order- 
of-magnitude change in step size shown in the figure. 
The central difference method uses the higher order 
operator and provides good accuracy over most of 
the step size range shown in the figure. The forward 
difference operator provides good accuracy with the 
smaller step sizes but begins to diverge earlier for the 
larger step sizes than the central difference method. 

Figure 5.29 shows the estimated derivative of dis- 
placement u\ at critical point 5 with respect to the 
rightmost span thickness /15 as a function of step size. 
In this case, the system mass matrix is a linear func- 
tion of this design variable and its derivative can be 
represented exactly by the forward difference oper- 
ator. The system stiffness matrix is a cubic func- 
tion of this design variable and its derivative can only 
be approximated by the forward difference operator. 
Still, the derivative approximation computed by the 
semianalytical method is very accurate except for the 
largest step size and is no worse for this case than the 
much more costly central difference method. Again, 
the forward difference operator results in substantial 
errors for the larger step sizes. 

Because this example has a relatively small num- 
ber of degrees of freedom, there is little condition 
error when small step sizes are used. To assess the 
effects of condition error which would occur for larger 
problems, the derivative approximations for the five- 
span beam problem were also calculated with 32- 
bit floating point precision compared with the 60- 
bit precision used in the studies described above. 
The estimated derivative of displacement, u 2 at criti- 
cal point 5 with respect to the point mass is plot- 
ted as a function of finite difference step size in 
figure 5.30. Derivative approximations are calculated 
using the semianalytical method, the central differ- 
ence method, and the forward difference method. For 
the larger step sizes, the results from all three meth- 
ods agree well with those calculated with 60-bit pre- 
cision. For step sizes smaller than 10 -4 in., there is 
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considerable error in the derivative approximations 
calculated with the forward and central difference 
methods. The derivative approximations calculated 
with the semianalytical method, however, are highly 
accurate over the entire range of step sizes shown in 
the figure. 



Figure 5.30. Effect of finite difference step size on accuracy 
of displacement derivative with respect to point 
mass design variable. Calculations performed 
with 32-bit precision; ramp load; undamped 
beam. 

5. 1.2. 3. Modal convergence of sensitivities. The 
first case considered is the undamped beam with the 
ramp load. Figure 5.31 shows the convergence of 
selected estimated derivatives of displacements with 
respect to m at various critical points. The mode 
displacement method was used and the derivative ap- 
proximations were calculated with the central differ- 
ence operator with updated modes. The convergence 
is good although slightly poorer than the conver- 
gence of the displacements themselves (fig. 5.8). The 
convergence of the estimated displacement deriva- 
tives with respect to the thickness design variable is 
similar. 

Although the modal convergence of the velocities 
for this case is good (fig. 5.9), the convergence of se- 
lected estimated derivatives of velocity with respect 
to m is generally poor (fig. 5.32). Given the poor 
convergence of the accelerations shown in fig. 5.10, 
it is not surprising that the convergence of the sen- 
sitivities of the accelerations is also very poor. From 
figure 5.33 it can be seen that the derivative approx- 
imations of the four selected critical point acceler- 
ations with respect to the thickness design variable 
are essentially not converging with increasing num- 
ber of modes. It should be pointed out again that 
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Figure 5.31. Modal convergence of derivatives of selected 
displacements with respect to mass design vari- 
able. Ramp load; undamped beam: mode dis- 
placement method; central difference operator. 

these derivative approximations of velocity and ac- 
celeration are calculated with the central difference 
method and updated mode shapes; thus, the numer- 
ical errors are minimized. The poor convergence ex- 
hibited in figures 5.32 and 5.33 is due to the poor ap- 
proximation of the sensitivities by the mode shapes. 

Similar modal convergence behavior is observed 
for sensitivities of the stress resultants. This is con- 
sistent with the poor convergence of the stress resul- 
tants calculated with the mode displacement method 
(figs. 5.11 and 5.12). Figure 5.34 shows the poor 
convergence of derivative approximations of selected 
bending moments with respect to the thickness de- 
sign variable. It can be seen that the convergence of 
the bending moment derivative approximation in the 
rightmost span with respect to the thickness in the 
rightmost span [dM^/dh^f) is especially poor. 

It was shown in the previous section that sev- 
eral approaches are available for overcoming the 
poor convergence of bending moments and shear 
forces in this beam example. The mode acceleration, 
static mode, and RWL methods all produced good 
modal convergence of bending moments and shears 
as shown in figures 5.14, 5.15, 5.16, and 5.18. Un- 
fortunately this type of dramatic improvement does 
not occur for the sensitivities of the stress resultants. 
Figure 5.35 shows the convergence of the bending 
moment derivative approximations with respect to 
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Figure 5.32. Modal convergence of derivative approxima- 
tions of selected velocities with respect to mass 
design variable. Ramp load; undamped beam; 
mode displacement method; central difference 
operator. 
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Figure 5.34. Modal convergence of derivative approxima- 
tions of selected bending moments with re- 
spect to thickness design variable. Ramp load; 
undamped beam; mode displacement method; 
central difference operator. 
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Figure 5.33. Modal convergence of derivative approxima- 
tions of selected accelerations with respect to 
thickness design variable. Ramp load; un- 
damped beam; mode displacement method; 
central difference operator. 
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Figure 5.35. Modal convergence of derivative approxima- 
tions of selected bending moments with respect 
to thickness design variable. Ramp load; un- 
damped beam; RWL method; central difference 
operator. 
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the thickness design variable where the analysis was 
done with the RWL method. As in the studies dis- 
cussed that used the mode displacement method, 
the sensitivities were calculated by the central dif- 
ference operator with the basis vectors updated for 
the perturbed design. Convergence of dM^/dh^ is 
somewhat improved compared with the mode dis- 
placement case. Other quantities show convergence 
similar to the mode displacement case; none of these 
convergence histories can be described as good. Con- 
vergence of the shear force derivative approximations 
with the RWL method is considerably worse than for 
the bending moments. 
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Figure 5.36. Modal convergence of derivative approxima- 
tions of selected bending moments with re- 
spect to mass design variable. Ramp load; un- 
damped beam; RWL method; semianalytical 
formulation. 

The semianalytical methods have also been used 
for calculating sensitivities of stress resultants. Fig- 
ure 5.36 shows the convergence of bending moment 
derivative approximations with respect to the mass 
design variable calculated with the fixed-mode, semi- 
analytical method and RWL vectors. The conver- 
gence is very similar, especially for larger numbers 
of modes, to that of the central difference method. 
The mode acceleration, semianalytical method, and 
the semianalytical method with approximate d$>/dx 
were also tried. Again, the modal convergence curves 
had the same jaggedness as for previous cases. 


Considering these difficulties with modal conver- 
gence for the ramp-loaded cases, especially poor 
convergence would be expected for the step-loaded 
case. For the ramp- loaded case, the convergence 
of displacement derivative approximations with re- 
spect to the mass design variable w r as reasonably 
good (fig. 5.31). For the step- loaded case, the modal 
convergence of the same displacement sensitivities is 
poor as shown in figure 5.37. The modal convergence 
of higher order sensitivities (velocities, accelerations, 
and stresses) is extremely poor. 

Adding damping slightly improved the modal con- 
vergence of the response quantities but did not com- 
pletely alleviate the convergence problem. The re- 
sult for sensitivities is similar. Figure 5.38 shows the 
convergence of velocity sensitivities for the discretely 
damped case. The convergence is slightly improved 
over the undamped case shown in figure 5.32, but the 
curves are still fairly jagged. Convergence of sensi- 
tivities for the case with modal damping is also very 
similar to that in the undamped case. 
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Figure 5.37. Modal convergence of derivative approxima- 
tions of selected displacements with respect to 
mass design variable. Step load; undamped 
beam; mode acceleration method; semianalyti- 
cal method. 

5. 1.2. 4- Mesh convergence of sensitivities. Just 
as for the response quantities, additional insight can 
be obtained by considering the convergence of their 
sensitivities with increasing number of elements per 
bay. The case of the stress resultants will be con- 
sidered, since they were shown to converge well with 
mesh refinement (fig. 5.26) but poorly as a function 


30 


o 

du j tdm 

tc ’ 
sec 

1.77 

□ 

diijldm 

2.18 

A 

d 14 3 tdm 

1.11 

• 

du 5 ldm 

.28 

■ 

dii^tdm 

2.29 



Figure 5.38. Modal convergence of derivative approxima- 
tions of selected velocities with respect to mass 
design variable. Ramp load; discretely damped 
beam; mode acceleration method; semianalyti- 
cal method. 

of number of vibration modes used in the analysis. 
Figure 5.39 shows the convergence of derivative ap- 
proximations of shear force with respect to the mass 
design variable as a function of number of elements 
per bay. Surprisingly, the convergence is extremely 
poor. 

5. 1.2. 5. Fixed versus updated modes in sensitivity 
calculations. As mentioned, the computational cost 
of updating the vibration modes for the perturbed 
analyses is substantial. The question of whether 
the modes from an initial design can be used in a 
finite-difference-based procedure to calculate sensi- 
tivities of the transient behavior has received consid- 
erable attention in the literature. In reference 35, it 
was shown that there is a substantial difference in 
the derivatives of aircraft flutter speeds when fixed 
modes are used rather than the updated modes. In 
reference 33, however, there was little difference in 
the derivatives of damping ratios for the five-span 
beam when either fixed or updated modes were used. 
This was investigated here with the same five-span, 
undamped beam under the step load. As shown in 
figure 5.37 where the derivative approximations were 
calculated with the semianalytical, mode accelera- 
tion method, convergence with respect to the number 
of modes is very slow. Figure 5.40 sho\frs the modal 
convergence of derivative approximations of selected 
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Figure 5.39. Convergence of derivative approximations of se- 
lected shear forces with respect to mass design 
variable as function of number of elements per 
span. Ramp load; undamped beam; mode dis- 
placement method; central difference operator. 

displacements with respect to h§ calculated with for- 
ward difference procedures. Results with both fixed 
and updated vibration modes are shown. Again, the 
convergence as a function of number of modes is poor. 
However, for all three derivative approximations, the 
results are nearly the same for both the fixed mode 
and updated mode cases. 

5.2. Composite Delta Wing Example 

The second example considered is an aircraft delta 
wing with laminated composite cover skins taken 
from reference 36 and described in detail in refer- 
ence 37. The finite element model of this structure 
is shown in figure 5.41. Since the wing is geomet- 
rically symmetric about the midplane through its 
thickness (X-Y plane), only the upper half of the 
wing is modeled, and boundary conditions enforcing 
antisymmetric motion are imposed on the joints ly- 
ing in the X-Y plane. The wing is also cantilevered 
at the root. The model contains a total of 88 joints 
with a total of 140 unconstrained degrees of freedom. 
The webs in the wing are made of titanium and are 
modeled with 70 shear panel finite elements along 
with rod elements through the thickness of the wing. 
The cover skins are made of a moderate- modulus 
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Figure 5.40. Modal convergence of derivative approxima- 
tions of selected displacements with respect to 
mass design variable calculated with both fixed 
and updated vibration modes. Step load; un- 
damped beam; forward difference methods. 


(E\ = 21 x 10 6 psi), graphite-epoxy material with 
0°, ±45°, and 90° lamina where the 0° material runs 
spanwise along the wing. These cover skins are mod- 
eled with membrane finite elements; thus, only the 
total thicknesses (and not the stacking sequence of 
plies) of each lamina are important. The structural 
mass is 6003 lb, but most of the wing mass is due to 
a fuel mass of 93 650 lb distributed over the joints. 
The spatial distribution of the load is the same as the 
static load from reference 37 and is roughly equiv- 
alent to a 144-psf pressure load on the wing skin. 
A step loading function was used as the time func- 
tion for all cases. The lowest 10 vibration frequen- 
cies for the wing are shown in table 5.3. Damping 
is accounted for by assuming 0.5 percent of critical 
damping for all modes. 


Table 5.3. Lowest 10 Vibration Frequencies 
for Delta Wing 


Mode 

Frequency, Hz 

1 

2.055 

2 

2.765 

3 

4.104 

4 

4.913 

5 

5.920 

6 

6.944 

7 

7.451 

8 

8.421 

9 

9.583 

10 

9.880 



Figure 5.41. Finite element model of composite delta wing. 


5.2 A. Wing Dynamic Response 

The character of the dynamic response of the 
delta wing is considerably different than that of the 
five-span beam. Shown in figure 5.42 is a time history 
of acceleration at the wingtip. Although 64 modes 
were included in the analysis, it is evident from 
figure 5.42 that only the low-frequency modes are 
being excited. The same is true for stresses as shown 
in the time history of figure 5.43. Shown in figure 5.43 
is r® which is a typical shear stress in a web. As can 
be seen, there is a small amount of higher frequency 
response superimposed on the predominant response 
frequency. However, the time history exhibits none of 
the high-frequency response present in the five-span 
beam. 

In contrast to the five-span beam example, the 
modal convergence of all response quantities consid- 
ered for the delta wing is quite good. Shown in fig- 
ures 5.44 and 5.45 are modal convergence plots for 
selected accelerations and stresses at critical points 
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Figure 5.43. Time history of shear stress in web in region 6 
for delta wing. 

calculated with the mode displacement method. A 
converged solution is reached with approximately 
20 or less modes for all the response quantities 
shown. Convergence is also good for response quan- 
tities when the mode acceleration or RWL methods 
are used instead of the mode displacement method. 
Shown in figure 5.46 is a convergence plot for the 
same stresses shown in figure 5.45 but calculated with 
RWL vectors. 

5.2.2. Sensitivities of Wing Dynamic Response 

5.2.2. 1. Design variables. The design variable 
definitions are the same as those in reference 37 and 
are shown in figure 5.47. As can be seen in figure 5.47, 
the skin is broken up into 16 regions. In each region 
there are three design variables— the thickness of the 
0° lamina, the thickness of the 90° lamina, and the 
thickness of the ±45° lamina. These design variables 
will be denoted t l B where i denotes the region of the 
wing skin, and 6 is either 0°, 90°, or 45° depending 
on the lamina orientation. Also shown in figure 5.47 
are the 12 design variables controlling the thickness 
of the webs. These will be denoted t l w where i denotes 
the particular web region. 

In calculating sensitivities of various response 
quantities, only a small subset of these design vari- 
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wing. 

ables were considered. Specifically, derivatives of se- 
lected displacement, velocity, acceleration, and stress 
quantities were calculated with respect to £ 90 > ^0 1 

4(1 ■ and l W- 

5. 2. 2. 2. Effect of finite difference step size. Com- 
pared with the five-span beam example, the system 
matrices for the delta wing are larger and have a more 
complicated connectivity. Since many of the signif- 
icant operations in the transient response analysis 
operate directly on these matrices, there is consider- 
able potential for accumulating round-off error. This 
round-off error along with the truncation error in the 
finite difference expressions is a concern in selecting 
a step size for a finite difference approximation to a 
derivative. 

A study was performed to consider the effect of 
step size in the forward difference and central differ- 
ence methods for the delta wing. Figure 5.48 shows 
derivative approximations of the wingtip acceleration 
at critical points with respect to selected thickness 
design variables as a function of the finite difference 
step size used. As seen in the figure, the step size was 
varied by factors of 10 from 10 " 7 to 10" 2 . The cen- 
tral difference method was not used with the 10 
step size because the backward perturbation from 
the nominal design would result in negative mem- 
ber thickness. One significant observation is that the 
acceptable step size range for the forward difference 
method is small — approximately 2 decades. A sec- 
ond observation is that the behavior of the central 
difference method as a function of step size is surpris- 
ingly good. It is expected that, for larger step sizes 
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Figure 5.45. Modal convergence of selected stresses for delta 
wing calculated with mode displacement 
method. 
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Figure 5.46. Modal convergence of selected stresses for delta 
wing calculated with RWL method. 




Figure 5.47. Design variable definitions for delta wing 
example. 
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Figure 5.48. Effect of finite difference step size on accu- 
racy of tip displacement derivative approxima- 
tions calculated with forward and central dif- 
ference methods with fixed modes. Delta wing 
example. 
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(1CT 3 ), the central difference method results in less 
error than the forward difference method. The unex- 
pected superior performance of the central difference 
method for the smaller step sizes is probably due to 
the symmetry of the difference operator. The round- 
off errors that occur with the positive and negative 
perturbations tend to cancel each other and thus pro- 
duce the better than expected accurate values for the 
sensitivities. 

The situation is similar for selected stress sensitiv- 
ities shown in figure 5.49. Most of the curves for the 
forward difference case have a small acceptable step 
size range. This is especially obvious for the deriva- 
tive daft/dt] $ where 1(T 5 is the apparent choice for 
step size. It should also be mentioned that these cal- 
culations were performed by using 64-bit arithmetic. 
In the five-span beam example, the effect of step size 
on displacement derivatives was not as severe even 
though, for one case, these were calculated with pre- 
dominantly 32-bit arithmetic (fig. 5.30). 
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Figure 5.49. Effect of finite difference step size on accuracy 
of stress derivative approximations calculated 
with forward and central difference methods 
with fixed modes. Delta wing example. 


The simple approach of selecting a single step size 
for use with all response quantities and all design 
variables was used here. This approach has the 
obvious advantage of simplicity but very questionable 


validity for the forward difference method and this 
delta wing example. From figure 5.49, there is 
significant error in da^/dt^ if greater than 10 ^ is 
used as the step size. However, if less than 10“° is 
used instead, du tip /dt}j 6 is in error. 

As noted the central difference method improves 
the range of acceptable step sizes but at the cost 
of an additional analysis for each design variable. 
Alternatively, the semianalytical method is partic- 
ularly attractive for this delta wing example. The 
stiffness matrices of the membrane and shear panel 
finite elements are linear functions of the thickness 
design variables. Thus, large values of the step size 
can be used to effectively eliminate the round-off 
errors in generating the derivative approximations 
of the stiffness and mass matrices required for the 
semianalytical method. 

5. 2. 2. 3. Modal convergence of sensitivities. Unlike 
the five-span beam example, the modal convergence 
of the displacement, velocity, and acceleration deriva- 
tives for the delta wing example is good. As an exam- 
ple, consider the reference case of acceleration sensi- 
tivities calculated with the central difference method 
with updated modes shown in figure 5.50. For all 
derivative approximations, convergence is achieved 
with 32 or less modes. Modal convergence for accel- 
eration sensitivities is equally good when the simple 
forward difference method with fixed modes is used 
as shown in figure 5.51. Figure 5.52 shows the conver- 
gence of acceleration sensitivities calculated with the 
semianalytical method with RWL vectors instead of 
vibration modes. Convergence is also good although 
slightly poorer than when modes are used. For ex- 
ample, approximately 40 RWL vectors are required 
for a converged value of dunp/dt^ compared with 
approximately 32 vibration mode shapes. 

The modal convergence of stress derivatives, how- 
ever, depends dramatically on whether fixed or 
updated modes are used in the calculation. The refer- 
ence case with the central difference operator uses up- 
dated modes, and as shown in figure 5.53, the modal 
convergence for all stress sensitivities is very good. 
Also the convergence of the stress sensitivities with 
the forward difference operator with updated modes 
as shown in figure 5.54 is very good with 24 or less 
modes yielding a converged solution. However, when 
the forward difference operator with fixed modes is 
used the modal convergence of the stress sensitivi- 
ties is very poor as shown in figure 5.55. For one 
derivative approximation, dcr^/dt^, the convergence 
is fairly good with approximately 24 modes yield- 
ing a converged solution. Especially poor conver- 
gence is observed for da^/dt^ (derivative of stress 
in the lamina with respect to its own thickness) 
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Figure 5.50. Modal convergence of tip acceleration sensitiv- 
ities for delta wing calculated with central dif- 
ference method. 
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Figure 5.51. Modal convergence of tip acceleration sensitiv- 
ities for delta wing calculated with forward dif- 
ference method with fixed modes. 


where approximately 100 modes are required for 
convergence. 

Using the semianalytical method with fixed 
modes does not improve the modal convergence of 
the stress sensitivities. Figure 5.56 shows the modal 
convergence of the same stress sensitivities as in the 
previous figures but calculated with the semianalyti- 
cal method with RWL vectors. The convergence be- 
havior for each derivative approximation here is very 
similar to that for the forward difference method with 
fixed modes. 

However, when the basis vectors are assumed 
to vary with the design variables and the modified 
modal method (see section 4.2.2) is used to approx- 
imate d$/drr, the results are significantly improved. 
Figure 5.57 shows the modal convergence of the same 
stress derivative approximations as shown in pre- 
vious figures. Here, the convergence is good with 
only around 24 modes required for convergence of 
the stress sensitivities. 

It was mentioned in chapter 4 that the modal 
method for approximating d&/dx produces no im- 
provement in the values of transient response sen- 
sitivities. This implies that including the modes 
in the modified modal method (see eq. (4.10)) 
may also not significantly improve the transient re- 
sponse sensitivities. This implication was tested 
by studying the modal convergence of the stress 
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Figure 5.52. Modal convergence of tip acceleration sensitiv- 
ities for delta wing calculated with semianalyt- 
ical method with RWL vectors. 
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Figure 5.53. Modal convergence of selected stress sensitivi- 
ties for delta wing calculated with central dif- 
ference method. 
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Figure 5.54. Modal convergence of selected stress sensitivi- 
ties for delta wing calculated with forward dif- 
ference method with updated modes. 
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Figure 5.55. Modal convergence of selected stress sensitivi- 
ties for delta wing calculated with forward dif- 
ference method with fixed modes. 
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Figure 5.56. Modal convergence of selected stress sensitivi- 
ties for delta wing calculated with semianalyti- 
cal method with RWL vectors. 
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Figure 5.57. Modal convergence of selected stress sensitiv- 
ities for delta wing calculated with semiana- 
lytical method with d&/dx approximated with 
modified modal method. 
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Modal convergence of selected stress sensitiv- 
ities for delta wing calculated with semiana- 
lytical method with d&/dx approximated with 
only pseudostatic solution. 


sensitivities with the use of the modified modal 
method but approximating d&/dx with only the 
pseudostatic term in equation (4.10). These results 
are shown in figure 5.58. Comparing this figure with 
figure 5.57 shows that for more than eight modes 
the results are nearly identical. It appears that a 
cheap, effective approximation to d&/dx in the semi- 
analytical formulation can be obtained with only the 
pseudostatic term from the modified modal method. 

For the five-span beam example, the convergence 
of the stresses was substantially improved by includ- 
ing the static solution via either the mode accelera- 
tion method, the static mode method, or the RWL 
method. The RWL method is attractive because it is 
cheaper to calculate n r RWL vectors than n r vibra- 
tion mode shapes. However, incorporating the modi- 
fied modal method in the sensitivity calculations with 
RWL vectors would seem to be impossible because 
it is derived to calculate the derivatives of vibration 
eigenvectors. (See eq. (4.10).) Regardless, it seems 
like a worthwhile numerical experiment to try using 
RWL vectors along with the pseudostatic correction 
term from the modified modal method in the variable 
mode, semianalytical formulation. One legitimate ar- 
gument for doing this is the well-known observation 


that the basis spanned by the RWL vectors is an ex- 
cellent approximation to the basis spanned by the 
eigenvectors. The results of this experiment for the 
modal convergence of the stresses in the delta wing 
are shown in figure 5.59. The convergence here is 
quite good also. For small numbers of modes the 
convergence is a bit erratic but in all cases the re- 
sults are good for more than 32 modes. The benefit 
of combining the RWL vectors with the pseudostatic 
approximation to d&/dx is that the RWL vectors 
add the often important static displacement compo- 
nent to the basis, whereas the pseudostatic term adds 
components reflecting the change in the design vari- 
able to the basis. 

As mentioned, the benefit of the mode accelera- 
tion method is that it also includes this pseudostatic 
term. The semianalytical, mode acceleration, sensi- 
tivity method described in chapter 4 was also applied 
to this delta wing example. Again, the modal conver- 
gence of the stress sensitivities shown in the previous 
figures is considered. Figure 5.60 shows the excel- 
lent convergence of the stress sensitivities. Clearly, 
the mode acceleration method provides the same im- 
provement in stress sensitivities as the semianalyti- 
cal method with a modified modal approximation to 
d&/dx. 
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Figure 5.59. Modal convergence of selected stress sensitivi- 
ties for delta wing calculated with seniianalyti- 
cal method with RWL vectors and pseudostatic 
approximation to d*P/dx. 
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Figure 5.60. Modal convergence of selected stress sensitivi- 
ties for delta wing calculated with semianalyti- 
cal mode acceleration method. . 


5.3. Stepped Cantilever Beam Example 

The third example considered is a cantilever beam 
with five different rectangular cross sections along the 
length. (See fig. 5.61.) This example is taken from 
reference 38 where its minimum mass design under 
a static tip load was considered. The thickness and 
width of the beam cross section in each of the five 
sections are given in the table insert on figure 5.61 
and represent an optimized design from reference 38. 
The beam is 200 in. long and, in the nominal case, 
each of the five sections has the same length. The 
material properties for the beam are also shown in 
figure 5.61. 
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Figure 5.61. Stepped cantilever beam with applied rota- 
tional acceleration at root. 


In most of the analyses, the beam is modeled with 
three finite elements per section. The transverse 
displacement and rotation are the nodal unknowns 
resulting in a total of 30 degrees of freedom for this 
case. The effect of different numbers of elements per 
section on the lowest 10 beam natural frequencies 
(with the beam clamped at the root) is shown in ta- 
ble 5.4. In the transient response analyses 0.5 per- 
cent of critical damping is included for each mode. 


5.3.1. Loading 

The loading for this stepped beam example is sig- 
nificantly different than for the first two examples. 
First, the load results from prescribing the accelera- 
tion at the beam root rather than by applying a force 
to the beam, and second, the time history as shown 
in figure 5.61 is more complicated than the simple 
step and ramp histories in the previous examples. 
The objective of this particular loading condition is 
to simulate the rotation of an appendage attached 
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Table 5.4. Lowest Frequencies for Stepped Cantilever Beam 


Mode 

Frequency, Hz, for 

3 elements 

4 elements 

5 elements 

6 elements 

1 

22.67 

22.67 

22.67 

22.67 

2 

102.67 

102.66 

102.66 

102.65 

3 

249.72 

249.62 

249.80 

249.55 

4 

440.57 

440.04 

439.80 

439.67 

5 

652.50 

650.82 

650.04 

649.62 

6 

878.48 

874.48 

872.58 

871.54 

7 

1093.36 

1086.15 

1082.61 

1080.64 

8 

1296.61 

1285.63 

1279.94 

1276.72 

9 

1479.81 

1465.74 

1457.74 

1453.09 

10 

1641.44 

1625.75 

1615.41 

1609.19 


to a relatively large mass (e.g., robotic arm). The 
acceleration history in figure 5.61 rotates the root of 
the beam through 10° in 0.18 sec. After 0.18 sec, 
the beam root is motionless while other points in the 
beam are undergoing dynamic motion. Beam dis- 
placements, velocities, and accelerations in the fol- 
lowing sections are with respect to the rotating coor- 
dinate system. 

This type of applied acceleration can be handled 
as an equivalent external force given as 

p = -Mrfl(i) (5.1) 

where r is a vector describing the linear rigid body 
rotation of the beam about its root and g(t) is the 
prescribed acceleration history given in figure 5.61. 
It should be noted that the applied force in this case 
depends on the system mass matrix; this must be 
considered in the sensitivity calculations. 

5.3.2. Stepped Beam Dynamic Response 

The transient behavior of the beam is strongly 
affected by the period of the loading. From table 5.4, 
the period of the lowest vibration mode is 0.044 sec, 
whereas the period of the square-wave loading is 
0.18 sec. From figure 5.62, it can be seen that in the 
time history of the beam tip displacement, this first 
mode predominates, and almost exactly four cycles 
occur during the period of the loading. After the 
load is removed, the displacement response at the 
tip is relatively small. The bending stress at the 
root has a time history similar to that of the tip 
displacement as can be seen in figure 5.63 but with 
slightly more participation from higher frequency 
modes. As expected, the acceleration time history 
for the tip as shown in figure 5.64 is considerably 
more jagged; this indicates the participation of many 


higher frequency modes. This behavior is largely due 
to the abrupt changes in loading in the square-wave 
input. Significant accelerations exist at the tip after 
the loading is removed. 



Figure 5.62. Time history of tip displacement for stepped 
cantilever beam. 
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Figure 5.63. Time history of root stress for stepped can- 
tilever beam. 



Figure 5.64. Time history of tip acceleration for stepped 
cantilever beam. 


40 



5. 3.2.1. Modal convergence. The first convergence 
study considered the effect of the number of finite 
elements per section on the convergence of the critical 
point displacements, velocities, and accelerations at 
the beam tip and stresses at the root. For all these 
quantities, the convergence is excellent. For example, 
the peak acceleration changes by less than 1 percent 
when the number of finite elements per section is 
varied from 3 to 8. Then, for the beam modeled 
with three elements per section, the effect of the 
number of modes used in the analysis was considered. 
Generally the convergence was better than expected. 
Figure 5.65 shows the modal convergence of the 
tip acceleration at two different critical time points 
calculated with the mode displacement method. The 
values are essentially converged with five modes. 
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Figure 5.65. Modal convergence of critical point tip acceler- 
ations for stepped cantilever beam. Mode dis- 
placement method. 


The modal convergence of the stress at the beam 
root is also rapid. Figure 5.66 shows the convergence 
of the root stress at two critical points calculated 
with the mode displacement method. No more than 
five modes are required for convergence. It was 
mentioned that there is a strong static component 
in the beam response during the period while the 
load is applied. Usually this requires the use of 
the mode acceleration method or RWL vectors for 
acceptable convergence of the stresses. Evidently the 
lowest vibration mode is close enough to the static 


displacement shape for this cantilever beam that the 
mode displacement method gives good values for the 
stresses. 

5 . 3. 2. 2. Use of RWL vectors in analysis. In 
the stepped beam and delta wing examples, the con- 
vergence with RWL vectors in analysis and sensitiv- 
ity calculations was generally as good or better than 
with vibration modes. The modal convergence in the 
stepped cantilever beam example when RWL vectors 
are used is very good also as seen in figure 5.67 for 
accelerations. 

As can be seen in figure 5.67, the largest number 
of RWL vectors used in the analysis is 20. In 
the convergence studies considering vibration modes 
(e.g., fig. 5.65), the full set of 30 modes was used. A 
complete set of RWL vectors could not be generated 
for this example because of ill-conditioning inherent 
in the numerical process (eqs. (2.26) through (2.29)). 
As additional RWL vectors are generated, round- 
off errors cause the vectors to become less and less 
orthogonal. Eventually, the vectors become linearly 
dependent; this results in a singular reduced system. 
In most practical applications of this R\\ L method, 
this singularity problem would not occur because the 
number of RWL vectors generated would be much 
smaller than the total number of degrees of freedom. 

5.3.3. Sensitivities of Stepped Beam Dynamic 

Response 

5.3.3. 1. Design variables. Two different classes of 
design variables are considered in this example. The 
first class is the set of beam thicknesses in each of 
the five sections. They are denoted hi, where i is the 
section number from figure 5.61. These are similar to 
thickness design variables considered in the five-span 
beam and delta wing examples. Sensitivity results 
are presented in the next sections with h\ and /15 
considered from this set. 

The second class of design variables is the set of 
lengths of the five sections in the beam. The beam 
length is fixed at 200 in.; thus, only four design vari- 
ables determine the lengths of the five sections. The 
four design variables are denoted /*, where l{ is the 
distance from the beam root to the end of the ith sec- 
tion. Sensitivity results are presented in the next sec- 
tions with l\ and I 4 considered from this set. In the 
structural optimization field this type of design vari- 
able is often referred to as a “shape” design variable 
and is studied separately from member thickness- 
type design variables. A recent study (ref. 39) con- 
sidered the calculation of static response sensitivi- 
ties with respect to shape design variables with the 
semianalytical method. It was found that numeri- 
cal difficulties in the semianalytical method resulted 
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Figure 5.66. Modal convergence of critical point root bend- 
ing stresses for stepped cantilever beam. Mode 
displacement method. 
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Figure 5.67. Modal convergence of critical point tip accel- 
erations for stepped cantilever beam. RWL 
vectors. 


in very large errors in sensitivities. This difficulty is 
addressed in sections 5.3.3.2 and 5. 3. 3. 3 for the tran- 
sient case. 
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Figure 5.68. Effect of finite difference step size on accu- 
racy of tip displacement derivatives with re- 
spect to thickness design variables for stepped 
cantilever beam. Forward and central differ- 
ence operators. 


5. 3. 3. 2. Effect of finite difference step size. 
It is shown in this section that, practically, the 
selection of finite difference step size is not a con- 
cern for this stepped beam example. A series of 
studies was performed to consider the effect of step 
size on both thickness and length sensitivities calcu- 
lated with finite difference and semi analytical meth- 
ods. The finite element model with three elements 
per section was used and all 30 modes were included. 
Figure 5.68 presents approximate derivatives of tip 
displacement with respect to section thicknesses cal- 
culated with overall forward and central difference 
methods. A key point to be made is that both meth- 
ods give excellent results for approximately an 8- 
decade step size range. For the large step size of 
10 1 in., the central difference operator generally 
gives better results than the forward difference oper- 
ator as would be expected. The results are nearly as 
good for sensitivities of the root stress with respect 
to the section thicknesses as shown in figure 5.69. 


42 


Compared with figure 5.68, there is slightly more er- 
ror for the smallest and largest step sizes but the 
sensitivities are still accurate over a very broad range 
of step sizes. If sensitivities of stresses with respect 
to the length design variables are considered, the re- 
sults are also very good. Figure 5.70 shows sensitivi- 
ties calculated with the forward and central difference 
methods. Again there is a broad range of step sizes 
that provide accurate sensitivities. For the smaller 
step sizes the results are generally less accurate than 
in figure 5.69, but for the 10~ 2 step sizes they arc 
more accurate. 
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Figure 5.69. Effect of finite difference step size on accu- 
racy of root stress derivatives with respect 
to thickness design variables for stepped can- 
tilever beam. Forward and central difference 
operators. 


It has been mentioned that severe numerical diffi- 
culties were uncovered in reference 39 when sensitiv- 
ities of static response were calculated with respect 
to shape design variables. The result of this numer- 
ical ill-conditioning could be seen by calculating the 
sensitivities with different finite difference step sizes 
used for approximating the derivatives of the stiff- 
ness matrix. For very small step sizes, the error in 
the sensitivities is due to round-off. For the larger 
step sizes, however, the errors in sensitivities were 
much larger than those due to truncation of the fi- 
nite difference operator and were found to be caused 


by basic ill-conditioning in solving the semianalytical 

equations. . . . 

This same phenomenon occurs when sensitivi- 
ties are calculated with a semianalytical method 
for the transient case. Figure 5.71 shows approxi- 
mate derivatives of root stress with respect to the 
length design variables calculated with the forward 
difference and semianalytical methods. Again, al 
30 modes are used in the analyses. For the smaller 
step sizes, the accuracy is significantly better for the 
semianalytical method than for the overall forward 
difference method. For the HP 2 step size, how- 
ever, the results from the forward difference method 
are excellent, whereas several of the sensitivities cal- 
culated with the semianalytical method exhibit ex- 
tremely large errors. This result is completely con- 
sistent with that in reference 39. Although in this 
example, there is a large range of step sizes where ac- 
curate sensitivities can be obtained, in general, tins 
would not be the case. Especially as the problem be- 
comes larger it is desirable to use larger step sizes in a 
semianalytical method, but this is severely restricted 
for shape design variables by the type of error shown 
in figure 5.71. 

5. 3. 3. 3. Modal convergence, of sensitivities. Most 
of the sensitivities exhibit the same good modal con- 
vergence as the response quantities. For example, 
the modal convergence of approximate derivatives of 
tip displacement with respect to hi and^/15 at differ- 
ent critical points is shown in figure 5.72. The sen- 
sitivities were calculated with the central difference 
method with updated modes and, as can be seen, the 
convergence is excellent. The convergence of tip ac- 
celeration derivative approximations is not as good 
as the displacement derivative approximations but is 
still acceptable as seen in figure 5.73. Again, these 
sensitivities are with respect to h\ and h b and are 
calculated with the central difference method with 
updated modes. 

Convergence is also good when sensitivities with 
respect to the length design variables are considered. 
Figure 5.74 shows the modal convergence of approxi- 
mate derivatives of acceleration with respect to h and 
; 4 calculated with the central difference method. A 
step size of 10 -5 was used to avoid the problem shown 
in figure 5.71. As can be seen in figure 5.74, conver- 
gence is achieved with approximately 10 modes. 

The modal convergence of stress sensitivities is 
similar to that for the delta wing example. When 
updated modes are used with an overall finite dif- 
ference method, the convergence is excellent. An 
example of this is shown in figure 5.75 where ap- 
proximate derivatives of the root stress with respect 
to hi and h b calculated with the forward difference 
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of approximate root stress derivatives with re- 
spect to length design variables for stepped can- 
tilever beam. Forward and central difference 
operators. 
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Figure 5.71. Effect of finite difference step size on accuracy 
of approximate root stress derivatives with re- 
spect to length design variables for stepped can- 
tilever beam. Overall forward difference and 
semianalytical methods. 
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Figure 5.72. Modal convergence of approximate derivatives 
of tip displacement with respect to thickness 
design variables for stepped cantilever beam. 
Mode displacement method; central difference 
operator. 
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Figure 5.73. Modal convergence of approximate derivatives 
of tip acceleration with respect to thickness 
design variables for stepped cantilever beam. 
Mode displacement method; central difference 
operator. 
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sec 

O dii ivp !dl x 0.09 
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Figure 5.74. Modal convergence of approximate derivatives 
of tip acceleration with respect to length de- 
sign variables for stepped cantilever beam. 
Mode displacement method; central difference 
operator. 

operator are shown. However, if fixed modes are used 
in a finite difference procedure, the modal conver- 
gence is much worse. Figure 5.76 shows an example 
of this for the same sensitivities as in figure 5.75. 
Also, if sensitivities of the root stress with respect 
to the length design variables (l\, I 4 ) are considered, 
the modal convergence is very poor. An example of 
this poor convergence is shown in figure 5.77. The 
convergence is similarly bad if the fixed-mode semi- 
analytical method is used instead of a finite difference 
method. Figure 5.78 shows the poor modal conver- 
gence of the same sensitivities as figure 5.77 but cal- 
culated with the fixed-mode, semianalytical method. 

For the delta wing example, remedies for the 
poor convergence of stress sensitivities in the semi- 
analytical method were based on approximating the 
mode shape derivatives d<i>/dx. These semianalyt- 
ical methods including approximations for d&/dx 
were also applied to this stepped beam exam- 
ple. First d&/dx was approximated with the 
modified modal method. The modal convergence 
of the stress sensitivities is now excellent as can 
be seen in figure 5.79. The degree of improve- 
ment can best be appreciated by comparing fig- 
ures 5.78 and 5.79 and noting that the range 
of the ordinate in figure 5.78 is much broader 
than in figure 5.79. Using only the first pseudo- 
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Figure 5.75. Modal convergence of approximate derivatives 
of root stress with respect to thickness design 
variables for stepped cantilever beam. Mode 
displacement method; forward difference oper- 
ator; updated modes. 

static term from the modified modal method as an 
approximation to d<f>/dx was also tried. As can 
be seen in figure 5.80, the convergence is adequate 
though not quite as good as when the complete mod- 
ified modal method is used. 

Just as in the delta wing example, a case was also 
considered where RWL vectors were used instead of 
vibration modes but their derivatives were computed 
by the modified modal method (version with pseudo- 
static term plus modes). Again, somewhat surpris- 
ingly, the modal convergence of the stress sensitivities 
is good as seen in figure 5.81. 

The semianalytical mode acceleration method 
was also tried as a remedy for the poor convergence 
of the stress sensitivities. Again, the very poor con- 
vergence is eliminated as can be seen in figure 5.82. 

5.4. Summary 

A number of different methods for calculating sen- 
sitivities of transient response quantities have been 
exercised on three example problems: a five-span 

beam, a composite aircraft wing, and a variable- 
cross-section beam. Two of the methods are over- 
all finite difference methods where the analysis is re- 
peated for perturbed designs. The other methods 
are termed semianalytical methods because they in- 
volve direct, analytical differentiation of the equa- 
tions of motion with finite difference approximations 
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Figure 5.76. Modal convergence of approximate derivatives 
of root stress with respect to thickness design 
variables for stepped cantilever beam example. 
Mode displacement method; forward difference 
operator; fixed modes. 
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Figure 5.78. Modal convergence of approximate derivatives 
of root stress with respect to length design 
variables for stepped cantilever beam. Mode 
displacement method; semianalytical method. 
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Figure 5.77. Modal convergence of approximate derivatives 
of root stress with respect to length design vari- 
ables for stepped cantilever beam. Mode dis- 
placement method; forward difference operator: 
fixed modes. 
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Figure 5.79. Modal convergence of approximate derivatives 
of root stress with respect to length design 
variables for stepped cantilever beam. Mode 
displacement method; semianalytical modified 
modal method. 
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Figure 5.80. Modal convergence of approximate derivatives 
of root stress with respect to length design 
variables for stepped cantilever beam. Mode 
displacement method; semianalytical one-term 
modified modal method. 


of the coefficient matrices. All the methods use basis 
vectors to reduce the dimensionality of the problem. 
Accordingly, the convergence of both the transient 
response quantities and their sensitivities as a func- 
tion of number of basis vectors was a key concern in 
this chapter. 

In the delta wing and stepped cantilever beam 
examples, the convergence of the response quanti- 
ties was consistently very good. However, this was 
not true with the five-span beam. With the five- 
span beam under a concentrated end moment and 
ramp time history, the convergence of displacements 
and velocities was adequate. However, the conver- 
gence of accelerations was poor. The convergence 
of stress resultants for this example depended on 
how they were calculated. When the mode displace- 
ment method was used, the convergence was quite 
poor. However, when the mode acceleration method, 
the Ritz-Wilson-Lanczos vector method, or the static 
mode method was used, the convergence was good. 
In cases where convergence was poor for the five-span 
beam, the addition of modal or discrete damping im- 
proved the convergence somewhat. However, it did 
not eliminate the convergence problems. 

The modal convergence of the sensitivities in the 
three examples is consistent with the convergence of 
the response quantities themselves. For the delta 
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Figure 5.81. Modal convergence of approximate derivatives 
of root stress with respect to length design vari- 
ables for stepped cantilever beam. RWL vec- 
tors; semianalytical modified modal method. 


O da TOO{/ ut x 


Jdl 

□ ^ a root^1 

A do ro Jdl x 

• do root ^4 

■ d<J TOOt fdU 


t c , 

sec 

0.07 

.11 

.16 

.02 

.16 



Figure 5.82. Modal convergence of approximate derivatives 
of root stress with respect to length design 
variables for stepped cantilever beam. Mode 
acceleration; semianalytical method. 
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wing and stepped cantilever beam examples, the 
convergence of sensitivities was generally good. For 
the five-span beam example, the convergence of dis- 
placement sensitivities was adequate but the conver- 
gence of velocities, accelerations, and stress resul- 
tants was generally poor. This poor convergence was 
observed for all the sensitivity calculation methods. 
Furthermore, it appears to be associated with the 
structure and loading because no improvement was 
observed as the number of finite elements per span 
was increased. 

In certain cases poor convergence of sensitivities 
was also observed for the delta wing and stepped can- 
tilever beam examples. When sensitivities of stresses 
were calculated with the fixed-mode overall finite 
difference methods or the fixed-mode semianalytical 
methods, the convergence was very poor. In large 
problems, however, updating the vibration modes in 
the overall finite difference methods or rigorously cal- 
culating derivatives of the mode shapes is very ex- 
pensive. The mode acceleration version of the semi- 
analytical method and the semianalytical method 
with mode shapes approximated with the modified 
modal method were devised to alleviate this poor 
convergence with lower computational cost. When 
both methods are applied to delta wing and stepped 
cantilever beam examples, the modal convergence of 
sensitivities is excellent. 

All the sensitivity calculation methods considered 
herein rely on finite difference operators. Thus step 
size selection is an important concern. The sys- 
tem stiffness and mass matrices are linear functions 
of many of the design variables in the three exam- 
ple problems. This allowed large step sizes to be 
used in the semianalytical methods to minimize the 
round-off errors and produce accurate derivatives of 
the stiffness and mass matrices. Also there is less 
opportunity for round-off error in calculating finite 
difference derivatives of just the coefficient matrices 
compared with finite difference derivatives of the 
overall response quantities. For these reasons, the 
semianalytical methods were consistently less sensi- 
tive to finite difference step size than the overall finite 
difference methods. 
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Chapter 6 

Computational Costs 

A consideration of the computational costs is es- 
sential for evaluating any numerical method. This is 
especially difficult in large-scale, finite-element-based 
procedures because there is often considerable “over- 
head” required in the practical implementation of a 
given numerical method. For example, most finite el- 
ement codes require only a small portion of a system 
matrix to be resident in central memory during fac- 
torization at any given time. The other portions of 
the matrix are read from disk and the factored por- 
tions are written to disk as required. A similar situa- 
tion can exist on virtual memory machines where the 
disk operations are transparent to the implementor. 
In these cases, the computer resources required are 
very implementation dependent. 

An approach that is common in the formal study 
of numerical methods is to evaluate the computa- 
tional cost by counting the number of floating point 
operations. There are some pitfalls to this approach. 
Sometimes, even for large problems, because of the 
required overhead it is impossible to achieve a prac- 
tical implementation that will execute as fast as the 
predictions from the operation count. At other times, 
especially on vector machines, it is possible for a 
method with a higher operation count to be faster 
than a method with a lower operation count. 

Nevertheless, this approach is used here, primar- 
ily to indicate the major trends in the costs of the 
methods, not to make fine distinctions between them. 
Following common practice, a floating point opera- 
tion (or “flop” ) is defined as the combination of a 
floating point multiply, add, and associated array in- 
dexing. In the rest of this chapter a floating point 
operation is often referred to simply as an operation. 

6.1. Costs of Basic Matrix Manipulations 

Multiplication of full matrices occurs in several 
places in the transient response and sensitivity meth- 
ods. The approximate number of floating point op- 
erations required to multiply a full l x m matrix and 
a m x n matrix is given as 

(6T) 


Solution of the reduced eigenproblem (eqs. (2.23)) 
is important in solving the system eigenproblem with 
subspace iteration (eqs. (2.3)) and in uncoupling the 
reduced system when basis vectors other than the 
eigenvectors are used. In both cases, it is necessary 
to solve a full, generalized eigenproblem for all n r 
eigenvalues and eigenvectors. Since eigenvalue solu- 
tion techniques are inherently iterative, the number 
of operations required for a converged solution can 
only be estimated. Reference 15 estimates the num- 
ber of operations for the complete solution of a gen- 
eralized eigenproblem with the Jacobi method as 

CVeig - 18 n* + 3671 (6.2) 

Other techniques for solving this eigenproblem may 
have a significantly different cost. However, it is 
shown that the cost of this eigensolution is small 
relative to other tasks in the sensitivity calculations. 

6.2. Costs of System Matrix Manipulations 

For the purpose of considering the computational 
costs of operations on system matrices (e.g., K, 
M), these matrices are considered to be stored in 
a banded form. In a banded form, only matrix ele- 
ments located near the diagonal are stored, the ma- 
trix elements outside this “bandwidth of the matrix 
are zero and are not stored or considered in oper- 
ations. Finite element problems yielding a stiffness 
matrix with a constant bandwidth are rarely encoun- 
tered in practice so most finite element codes use 
more sophisticated and efficient schemes for storing 
system matrices. However, having a single, easily 
understood number to characterize the sparsity of 
a system matrix (the bandwidth) is convenient in 
approximating computational costs. Although few 
finite element problems have precisely a constant 
bandwidth, this assumption is accurate enough in 
many cases to get reasonable estimates for a relative 
number of operations in a numerical procedure. 

From reference 40 the cost in number of opera- 
tions of factoring a banded system matrix of order 
n g is given as 

C bfac = \&{f* + 3)n s - y - ? ~ f 0 M 

where (3 is one half the bandwidth (excluding the 
diagonal). Also from reference 40, the cost of a 
single solution of a banded system, given the factored 
matrix, is given as 

Cbsol = 2 (£ + l)«p - + : ) ( 6 - 4 ) 


Cfmul = l mn 
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The cost of multiplying a banded system matrix and 
a single vector is given as 

0 )mul = (2/? + (6-5) 


6.3. Cost of Basis Reduction 

The process of reducing the degrees of freedom 
from rig to n r requires the matrix triple product op- 
erations shown in equations (2.6), (2.7), and (2.8). 
Since this process is used in all the sensitivity meth- 
ods, the cost is considered separately here. In per- 
forming this operation, n r system vectors are multi- 
plied by a system matrix. Then n r (n r + l)/2 inner 
products (for a symmetric system matrix) between 
system vectors are performed. The total number of 
required floating point operations is 

^red — n rCbmul H (6.6) 


6.4. Cost of System Eigensolution 

The cost of solving the generalized vibration 
eigenvalue problem is even more difficult to estimate 
than the cost of solving the reduced eigenproblem. 
The numerical techniques vary widely among differ- 
ent analysis codes. Furthermore, a technique used 
for one problem might be totally inappropriate when 
applied to a different problem. Nevertheless, some 
assumptions are made here that will hopefully lead 
to a reasonable estimate of computational costs for a 
fairly broad class of problems. 

First, it is assumed that the eigenvalue problem 
is solved with a subspace iteration technique with 
shifts (for example, ref. 15). In recent years, software 
based on this approach has become common. Also, 
the eigensolver, E4, in the EAL software used in this 
study (ref. 23) is based on this approach. It is also 
necessary to make assumptions about the number of 
vectors used in the subspace and the number of itera- 
tions at a given shift point required to converge some 
subset of these vectors to eigenvectors. The follow- 
ing numbers were used for these quantities with the 
realization that they may be optimum for only a few 
problems. Also, it is assumed that the eigensolution 
is being performed for a slightly perturbed model and 
the eigenvectors from the initial model are available 
as the initial subspace. At each shift point /it S s = 16 
vectors are included in the subspace. After n it = 2 
iterations, n css = 8 of these vectors have converged 
to eigenvectors. 


The number of shifts or number of factorizations 
required is approximately 

^shift — H 1 (6-7) 

™CSS 

At each iteration, the inverse power operation re- 
quires a matrix product between M and nt ss vec- 
tors followed by s solutions of the system equations 
based on the current factored K. The basis reduction 
operation for both K and M requires n tS s(^tss + 1) 
system vector inner products. Next, the eigen- 
problem of reduced order n tss must be solved. This 
cost is given in equation (6.2) with n r replaced by 
ntss- Finally, the updated set of approximate system 
eigenvectors must be formed as a linear combination 
of the current approximation. This requires n^ ss rig 
operations. The approximate cost of solving the^ys- 
tem eigenproblem for n r modes and frequencies can 
be written as 

Qig = n shiftQ)fac + n shift n it[rctssCbmul 

T ^tssCbsol T fltss(^tss + 1 )Tlg 

+ l^ n tss + 36n, ss + n ( 2 SK rig] (6.8) 


6.5. Cost of Generating RWL Vectors 

As has been demonstrated, RWL vectors are an 
attractive alternative to vibration mode shapes for 
basis reduction in transient response analysis. It 
has been mentioned previously that generation of the 
RWL vectors is considerably cheaper than vibration 
modes. An estimate of this cost in number of floating 
point operations is derived here. 

First, a factorization of the system K is required. 
The system equations are solved n r times based on 
the factored K. The generation of right-hand-side 
vectors requires n r — 1 matrix products between M 
and a vector. Another key step in the process is 
the Gram-Schmidt orthogonalization as indicated in 
equation (2.27). For all vectors, this requires ;z r - 1 
multiplications of a vector by M and n r (n r - l)/2 
vector inner products. The scaling of each vector 
requires n r vector inner products and n r divisions of 
a system vector by a scalar. Writing the total number 
of floating point operations in expanded form yields 

Crwl = c bfac + n r C bsol + 2(n r - l)C bmul 

n r (n r - 1) 

+ 5 rig + 2 n r n g (6.9) 
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6.6. Cost of Model Generation 

The generation of the finite element model re- 
quires processing of the input, forming elemental ma- 
trices, and forming global system matrices. Most 
of the sensitivity calculation methods require gen- 
eration of a single perturbed model for each design 
variable. The central difference method, however, re- 
quires the generation of two perturbed models. Thus 
to compare the central difference method with the 
other methods, an estimate of the model generation 
cost is required. This cost is difficult to calculate 
in general. For the purposes herein this cost is es- 
timated empirically with EAL by observing the exe- 
cution time for model generation relative to matrix 
multiplication for a number of models. From these 
experiments it was observed that the predominant 
element type in the model substantially affects the 
cost. That is, forming the element matrices in a 
model composed of three-dimensional solid elements 
is much more costly than in a model composed of rod 
elements. The estimate for model generation cost 

used here, . 

Cmodel = 100 Png (6.10) 

roughly approximates the cost for a model with 
two-dimensional, plate-tvpe elements in EAL but 
would be significantly in error for predominantly one- 
dimensional or three-dimensional models. 

6.7. Cost of Integration of Reduced System 


The basic operation for integrating the reduced 
system is shown in equation (2.30). The two ma 
trix multiplications shown in equation (2.30) are per- 
formed at every time step. If equations (2.5) are cou- 
pled, W,, and N (J are full and the explicit matrix 
multiplication must be performed. In this case, the 
total number of floating point operations for integra- 
tion of the system is given as 

Cinte = 8 nlm (6.11) 


where n* is the number of time steps in the analy 
sis If equations (2.5) are not coupled, W jj and Njj 
are diagonal, and this fact can be exploited to sub- 
stantially reduce the cost of integrating the system. 
The number of floating point operations in this case 


is given as 


Cinte = 8 n r n t 


( 6 . 12 ) 


When the number of equations in the reduced sys- 
tem n r is large, the difference between C inte in 
equations (6.11) and (6.12) is very large. For the 
comparisons of sensitivity methods in this chapter, 
equation (6.12) is used to estimate the integration 


cost. When vectors other than vibration modes are 
used or vibration modes for an initial model are used 
with a perturbed model, the equations are first un- 
coupled by solving the reduced-order eigenproblem. 

6.8. Cost of Back Transformation for 
Physical Response Quantities 

After the reduced equations have been solved, it 
is necessary to recover the physical displacements, 
velocities, accelerations, and stresses (or stress resul- 
tants) of interest. Usually the quantities of interest 
are only a subset of all possible quantities available 
from the finite element model. In the critical point 
constraint formulation described in chapter 3 it is 
necessary to recover the physical response quantities 
only at the critical times. That is, the back trans- 
formation is performed at only 5 to 20 critical points 
rather than at thousands of time steps. The cost of 
the basic back transformation operation is 

E’back = npUrUc (6-13) 

where n„ is the number of physical quantities being 
recovered from the modal values and n r is the number 
of critical points. The costs of back transformation 
in the specific sensitivity methods will be expressed 
as a multiple of this basic cost. 

6.9. Cost of Sensitivity Calculation 

Methods 

Because all the sensitivity calculation methods 
require the dynamic analysis of the initial model, this 
component of the cost can be neglected in comparing 
the different methods. In addition, in all the methods 
the basic operations are repeated for each design 
variable so the costs estimated below are per design 
variable. Also, to simplify the cost analysis, the 
models are assumed to be undamped so that any 
operations dealing with modal damping or system 
damping matrices are not included. 

6.9.1 . Finite Difference Methods 

Both forward and central difference methods for 
calculating sensitivities were considered in chapter 4. 
In the central difference method, the basic operations 
of the forward difference method are performed twice; 
therefore the cost is approximately twice that of the 
forward difference method. Costs are derived here 
for the forward difference method. In both finite 
difference methods, the basis vectors can be the 
same as for the original model (fixed) or recalculated 
for the perturbed model (updated). The cost with 
updated modes presented herein is based on using 
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natural vibration modes. An alternative of using 
RWL vectors is considered separately. 

The first step in the forward difference method is 
evaluation of the perturbed model. If the modes are 
being updated, the eigenproblem is solved. Other- 
wise, the original modes are used to reduce the basis, 
and the reduced-order eigenproblem is solved to un- 
couple the transient equations. The uncoupled equa- 
tions for the perturbed system are then integrated 
and the n p physical quantities calculated at the n c 
critical time points. The cost of the actual difference 
operation is very small and is therefore neglected. 
For the fixed-mode case, the total cost is 

Cfdfix = C'model + 2C re( j + C reig + riyUg 

+ Cinte + Cback (6. 14) 

For the updated mode case, the total cost is 

Qdupd = Crnodel + Ceig + Cj u te + C’) )ac k (6.15) 

6.9.2. Semianalytical Method With Fixed Modes 

The semianalytical method begins by evaluating 
the perturbed model. Then dM/dx and dK/dx are 
formed using a forward difference operator. Each 
derivative requires about /3n g operations. Then 
the basis reduction operation is applied to both 
derivative matrices. Formation of the right-hand- 
side pseudo load (eqs. (4.6)) is a fairly costly op- 
eration and the two matrix products (dM/dx)q and 
(dK/dx)q require about n 2 n t operations each. Fi- 
nally, the uncoupled equations are integrated and 
the physical sensitivities recovered. For the purposes 
of cost estimation, a single quantity n p is used as 
the total number of required physical sensitivities. 

In the semianalytical methods, however, the specific 
procedure for recovering the sensitivities depends on 
whether the quantity is a displacement, velocity, ac- 
celeration, or stress sensitivity. In estimating the 
costs of this back transformation operation these dif- 
ferences are ignored. One justification for this ap- 
proach is that the cost of back transformation is usu- 
ally small relative to other costs in the sensitivity 
calculation. In this fixed-mode, semianalytical 
method, approximately the same number of opera- 
tions is required for the recovery of physical sensitiv- 
ities as in the finite difference methods. The total 
number of floating point operations can be written 
as 

^safix = Cmodel + ‘2 3 Tig + 2 C rec ( + 2 n 2 tlt 

+ Qnte + Chack ( 6 . 16 ) 


6.9.3. Semianalytical Method With Approximate 
d<F/dx 

Just as in the fixed-mode, semianalytical method, 
evaluating the perturbed model and forming dM/dx 
and dK/dx is the first step. The next step is using 
the modified modal method to approximate d<F /dx. 

The procedure for the modified modal method is 
given in equations (4.9), (4.10), (4.11), and (4.12). 
The calculation of the n r pseudostatic contributions 
requires the formation of n r right-hand-side vectors 
and n r solutions of the system equations. The forma- 
tion of the A jk participation factors requires approx- 
imately n r system matrix additions plus the equiv- 
alent of a triple product basis reduction operation. 
Forming the linear combination of pseudostatic term 
and eigenvectors requires n r n g operations. The total 
cost in number of floating point operations for the 
modified modal method is 

^mmod = 2 n r n g -)- n r C| )SO ] -)- n r 0n g 

+ C Ted + n r n g (6.17) 

Given d<fr/dx, the derivatives of the reduced sys- 
tem matrices can be formed. For both dM/dx and 
dK/dx, two triple product, basis reductions plus 
n r vector inner products (for the d$ r /dxM$ term 
since is already available) are required as shown 
in equation (4.8). The right-hand-side formation and 
integration of the reduced sensitivity equations are 
identical to the fixed-mode semianalytical method. 
Because of the nonzero d<&/ dx, recovery of the physi- 
cal sensitivities is more complicated than in the fixed- 
mode case. Approximately twice the number of op- 
erations is required in the back transformation since 
both $ and d<F/ dx terms must be considered as 
shown in equations (4.13). The total cost for the 
variable-mode semianalytical method is 

Qaiipd = ^ ’model + 2f3n g + C mmo( i -f 4C re( j + 2 n r n g 

+ 2n'j.n t + C inte + 2Q, ack (6.18) 


6.9.4. Semianalytical Mode Acceleration Method 

Since dq/dx and dq/dx are obtained from the 
fixed-mode semianalytical method, the operations in 
equation (6.16) (except C back ) are required in apply- 
the mode acceleration method. The back trans- 
formation operations for displacement and stress 
sensitivities are more complicated as seen in equa- 
tions (4.19) and (4.20). The cost of forming the co- 
efficients in equations (4.19) is dominated by multi- 
plying a vector by a system matrix, adding n r + 1 
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system vectors, and solving the system equations for 
n T + 1 pseudostatic vectors. Again, the assumption 
is made that the model is undamped so that the C 
and the dC/dx terms in equations (4.19) are zero. 

The back transformation procedure for displace- 
ment and stress sensitivities involves application of 
equations (4.19) and (4.20) for each quantity at 
the critical times. Velocity and acceleration terms 
are calculated as in the fixed mode, semianalytical 
method. Again, with only a single quantity for the 
number of back transformed quantities rip, the cost 
can only be roughly estimated as 4Cback- The to- 
tal cost for the semianalytical, mode acceleration 
method can then then be written as 

Csamacc = Cm 0c iel "b 2pTlg + 2C , rec j + 2u r 7lt 

+ Cbmiil + ( n r + + ( nr + ^^bsol 

+ Ci n te + ^Cback (6.19) 

6.10. Analysis of Cost For Various Models 

With the expressions for computational cost in 
the previous sections, it is now possible to evaluate 
the use of the sensitivity calculation methods on var- 
ious examples. The first three examples are those 
considered in chapter 5. These three examples, how- 
ever, are all rather small compared with the class 
of problems envisioned for the production use of the 
sensitivity methods. Accordingly, two other hypo- 
thetical problems with a larger number of degrees of 
freedom have been included. 

The key parameters from the five problems re- 
quired for the cost analysis are shown in table 6.1. 
Several points should be made about these parame- 
ters. The two beam problems have a small number 
of degrees of freedom and a very small bandwidth 
and, as a result, a small cost for system matrix fac- 
torization. This is unusual in finite element anal- 
ysis. Medium model A and large model B repre- 
sent a typical medium size linear dynamics problem 
and a rather large ambitious problem, respectively. 
Medium model A also is complicated by the fact 
that 100 vectors are assumed to be required in the 
transient analysis. In all five examples, a relatively 
large number of time steps are used in the transient 
analysis. 

6.10.1 . Cost of Computational Subtasks 

Table 6.2 shows the number of floating point op- 
erations required for different computational tasks 
for the five example problems. Examining the costs 


for these subtasks gives some clues to the costs of dif- 
ferent sensitivity calculation methods. For the first 
three examples, the cost of system matrix factoriza- 
tion is low. For the two larger hypothetical examples 
the factorization cost is much higher relative to that 
of other tasks. In the first three examples, the cost of 
integrating the reduced equations is substantial even 
though the equations are uncoupled. For models A 
and B, the integration cost is 1 to 2 orders of mag- 
nitude less than the other subtask costs in table 6.2. 
Consistently, in all five examples, the cost of perform- 
ing the triple product basis reduction is high. For the 
three small problems, this cost is significantly higher 
than the factorization cost. For medium model A, 
this cost is also much higher than the factorization 
cost, but this is primarily due to the requirement of 
100 vectors in the reduced system. Even in model B, 
however, the basis reduction cost is only a little less 
than one half the factorization cost. One conclusion 
is that the number of vectors in the reduced system 
substantially affects the cost of the analysis even if 
the vectors are not updated for the current model. 


Table 6.1. Parameters Governing Computational Costs 


Model 

n g 

0 

n r 

n* 

n P 

Ti c 

Five-span beam 

32 

3 

18 

6000 

25 

10 

Delta wing 

264 

30 

20 

30 000 

13 

5 

Stepped beam 

32 

3 

20 

30000 

4 

5 

Medium model A 

3000 

100 

100 

10000 

50 

10 

Large model B 

12 000 

300 

30 

20000 

200 

10 


The use of RWL vectors in the transient and sen- 
sitivity analyses was considered in chapter 5. Here, 
the cost of generating RWL vectors compared with 
vibration modes is considered. Table 6.2 shows the 
cost of system matrix eigensolution C e i g and RWL 
vector generation Cr\vl f° r ^ ve example prob- 
lems. In every case the generation of RWL vectors is 
cheaper than the eigensolution. In the beam exam- 
ples, Crwl is more than an order of magnitude less 
than Ceig. This results from the unusual situation in 
which the number of required eigenvectors is nearly 
the same as the total number of degrees of freedom. 
In this case, the solution of the reduced eigenproblem 
artificially raises the cost of the system eigensolution. 
The other three examples show C e i g to be three or 
four times Crwl- This is probably a much more ac- 
curate estimate of the cost savings obtained by using 
RWL vectors instead of eigenvectors. 
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Table 6.2. Number of Operations for Selected 
Computational Subtasks 


Model 

Cbfac 

b^red 

^eig 

Crwl 

Cjntc 

Five-span beam 

2.7 

X 

10 2 

9.5 

X 

10 3 

6.6 

X 

10 5 

1.8 

X 

10 4 

8.6 

X 

10 5 

Delta wing 

1.2 

X 

10 5 

3.8 

X 

10 5 

4.8 

X 

10® 

1.1 

X 

10® 

4.8 

X 

10® 

Stepped beam 

2.7 

X 

10 2 

1.1 

X 

10 4 

6.6 

X 

10 5 

2.1 

X 

10 4 

4.8 

X 

10® 

Medium model A 

1.5 

X 

10 7 

7.5 

X 

10 7 

7-4 

X 

10 8 

2.1 

X 

10 8 

8.0 

X 

10® 

Large model B 

5.4 

X 

10 8 

2.2 

X 

10 8 

4.0 

X 

10 9 

1.2 

X 

10 9 

4.8 

X 

10® 


6 A 0.2. Comparison of Costs for Five Sensitivity 
Methods 

The primary objective of this chapter is to com- 
pare the costs in number of floating point operations 
of the sensitivity methods. This is summarized for 
five sensitivity methods, for the five examples in ta- 
ble 6.3. It is believed that these five sensitivity meth- 
ods are all practical alternatives for large-order prob- 
lems. This belief is substantiated by the fact that for 
all five examples the difference among the five costs 
is less than 1 order of magnitude. 

Table 6.3. Overall Operation Costs for Five Sensitivity Methods 


Model 

b'fdfix 

^Tdupd 

c 

safix 

^saupd 

t^samacc 

Five-span beam 

1.0 

X 

10® 

1.5 

X 

10® 

4.8 

X 

10® 

4.8 

X 

10® 

4.8 

X 

10® 

Delta wing 

6.5 

X 

10® 

1.0 

X 

10 7 

3.0 

X 

10 7 

3.2 

X 

10 7 

3.1 

X 

10 7 

Stepped beam 

5.0 

X 

10® 

5.5 

X 

10® 

2.9 

X 

10 7 

2.9 

X 

10 7 

2.9 

X 

10 7 

Medium model A 

2.4 

X 

10 8 

7.8 

X 

10 8 

3.9 

X 

10 8 

6.8 

X 

10 8 

4.5 

X 

10 8 

Large model B 

8.2 

X 

10 8 

4.4 

X 

10 9 

8.5 

X 

10 8 

1.7 

X 

10 9 

1.1 

X 

10 9 


The forward difference method with fixed modes 
is consistently the cheapest method. However, this 
low computational cost must be weighted against 
the pitfalls of the method discussed in chapter 5. 
The cost of a fixed-mode central difference method 
which is approximately twice the forward difference 
cost would also be quite competitive with the other 
methods and would lessen the sensitivity to finite 
difference step size. For the two larger problems, 
the forward difference method with updated modes 
is relatively expensive and an updated-mode central 
difference method would be extremely expensive for 
larger problems. 

In the three smaller problems the semianalytical 
methods require significantly more operations than 
the finite difference methods. This is primarily be- 
cause the larger number of time steps makes the 
calculation of the right-hand-side pseudo load rela- 
tively large. For the two larger problems, however, 


the fixed- mode semianalytical method is quite com- 
petitive with the forward difference method. For 
model A, it is less than twice the cost of the finite 
difference method, and for model B, it is essentially 
the same. 

The number of basis vectors used is a key pa- 
rameter in both the analysis and sensitivity calcu- 
lations. Table 6.1 shows the number of modes used 
in the baseline cost analyses for the five examples. 
Here, the effect of the number of modes on the over- 
all sensitivity costs is considered. First, the delta 
wing example, which is representative of a typical 
small problem, is considered. Shown in figure 6.1 is 
the cost for the five methods plotted as a function 
of number of modes used. The number of modes 
ranges from 20 to 100. The values of the other pa- 
rameters in the problem are in table 6.1. The key 
result from figure 6.1 is that the semianalytical meth- 
ods are much more costly than the finite difference 
methods for large numbers of modes. There are two 
reasons for this: first, because the problem is small, 
calculation of the vibration modes is relatively cheap, 
and second, because there are a large number of time 
steps, formation of the right-hand side in the sensitiv- 
ity equations for the semianalytical methods is quite 
costly when the number of modes used is large. 

For the large model B example, the result of vary- 
ing the number of modes is very different. For this 
example, the cost of the five sensit ivity methods plot- 
ted as a function of number of modes is shown in fig- 
ure 6.2. In this example, the calculation of the modes 
is a very costly operation. Accordingly, the forward 
difference method with updated modes is substan- 
tially more costly than the other methods for large 
numbers of modes. The fixed-mode forward differ- 
ence and semianalytical methods show only moderate 
increases in cost as the number of modes is increased. 

It was mentioned above that a relatively large 
number of time steps are used in the five examples. 
The effect of the number of time steps on the overall 
sensitivity calculation costs is considered here. The 
delta wing example is considered as representative 
of a small problem; the large model B example, of a 
large problem. For the delta wing, the computational 
costs for the five sensitivity methods are plotted as 
a function of the number of time steps in figure 6.3. 
The values of the sensitivity calculation costs here 
are similar to those in figure 6.1; the forward differ- 
ence methods show only moderate cost increases for 
larger numbers of time steps and the semianalytical 
methods show substantial cost increases. Again, the 
reason is that the right-hand-side formation in the 
semianalytical methods is a substantial part of the 
total cost in small problems. 
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Figure 6.1. Cost of sensitivity calculation methods as func- 
tion of number of modes for delta wing example. 
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Figure 6.2. Cost of sensitivity calculation methods as func- 
tion of number of modes for model B example. 
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Figure 6.3. Cost of sensitivity calculation methods as func- 
tion of number of time steps for delta wing 


example. 
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Figure 6.4. Cost of sensitivity calculation methods as func- 
tion of number of time steps for model B 


example. 
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The cost results from a large problem, the 
model B example, are very different, however. The 
cost as a function of number of time steps is plot- 
ted in figure 6.4. Obviously, there is practically no 
change in cost for any of the methods as a function of 
number of time steps. For this large problem, both 
the cost of integrating the uncoupled equations and 
the cost of forming the right-hand side in the semi- 
analytical methods are 2 to 3 orders of magnitude 
less than the total cost. 

6.11. Summary 

The main objective of this chapter is a compari- 
son of the computational costs in number of floating 
point operations of the sensitivity calculation meth- 
ods. Five example problems were considered -the 
three example problems from chapter 5, which are all 
fairly small and two larger hypothetical examples. 

Many of the results depend significantly on 
whether the problem is one of the three smaller ex- 
amples or one of the two larger hypothetical exam- 
ples. In the three smaller examples, the cost of sys- 
tem matrix factorization is low, whereas in the larger 
problems, this cost is quite high. When the cost of 
factorization is high, the system eigenproblem is es- 
pecially costly. In the smaller problems, operations 
repeated for the reduced problem at each time step 
(such as integration of the uncoupled equations) are 
a significant percentage of the total sensitivity calcu- 
lation cost. For large problems, the relative cost of 
these operations is small. 

For all five examples, the forward difference 
method with fixed modes was the cheapest. For 
the smaller problems the forward difference method 
with updated modes had a relatively low cost, but 
for the larger problems the cost was quite high. For 
the larger problems the semianalytical method with 
fixed modes and the semianalytical mode accelera- 
tion method have costs that are relatively competi- 
tive with the fixed-inode forward difference method. 

In all cases, the semianalytical method with approx- 
imate eigenvector derivatives was one of the more 
costly methods. 

It was shown in chapter 5 that for two exam- 
ples the accuracy of the stress sensitivities for small 
numbers of basis vectors was extremely poor. It was 
demonstrated that the semianalytical mode accelera- 
tion method was one means of dramatically improv- 
ing this accuracy. From the results of this chapter, 
the semianalytical mode acceleration method is only 
slightly more costly than the fixed-mode forward dif- 
ference and semianalytical methods. Given the un- 
acceptable accuracy of these fixed-mode methods for 


two of the examples, the semianalytical mode accel 
eration method appears to be the method of choice 
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Chapter 7 

Concluding Remarks 

Several methods have been developed and eval- 
uated for calculating sensitivities of displacements, 
velocities, accelerations, and stresses in linear, struc- 
tural, transient response problems. Two of the meth- 
ods are overall finite difference methods where the 
analysis is repeated for perturbed designs. The other 
methods are termed semianalytical methods because 
they involve direct analytical differentiation of the 
equations of motion with finite difference approxima- 
tions of the coefficient matrices. The different sen- 
sitivity methods were evaluated by applying them 
to three example problems: a five-span simply sup- 
ported beam loaded with an end moment, an aircraft 
wing loaded with a distributed pressure, and a can- 
tilever beam with a stepped cross section loaded with 
an applied root angular acceleration. 

An important issue in calculating transient re- 
sponse sensitivities for use in formal optimization 
procedures is how to define the constraints. Two 
common approaches are to integrate the response 
quantity over time or to pick the maximum (or min- 
imum) value of the response quantity in time. Both 
these approaches have drawbacks. An alternative 
critical point constraint approach was implemented 
which identifies the most important response points 
along the time history. A method for identifying 
these critical points was devised that, based on the 
three examples considered, appears to be very effec- 
tive even for very jagged response histories. 

All the analyses and sensitivity methods consid- 
ered use approximation vectors to reduce the number 
of degrees of freedom in the analysis. Vibration mode 
shapes, Ritz-Wilson-Lanczos vectors, and static dis- 
placement shapes were used in the analysis and sen- 
sitivity calculations. The key question when an ap- 
proximate reduced basis is used in an analysis is how 
many basis vectors are required for an accurate ap- 
proximation to the finite element solution. It was 
generally found that, if the accuracy of the response 
quantities was poor, the accuracy of the sensitivities 
was extremely poor. In a number of cases, however, 
even though the accuracy of the response quantities 
was adequate, the accuracy of sensitivities was poor. 
This is discussed further below. In all cases consid- 


ered herein, the accuracy as a function of the num- 
ber of vectors for both the response quantities and 
sensitivities with Ritz-Wilson-Lanczos vectors was as 
good or better than with vibration modes. Since the 
generation of Ritz-Wilson-Lanczos vectors is cheaper 
than vibration modes, they appear to provide a more 
cost-effective alternative to modes in many cases. 

A goal in considering sensitivity methods in this 
study is that they be suitable for very large-order fi- 
nite element analysis. In these types of problems, 
a complete vibration analysis tor each perturbed 
model is impractical because of the high computa- 
tional cost. To reduce this cost, one approach which 
was studied herein is to use the basis vectors from 
the initial model to approximate the response in the 
perturbed model. This often provides an effective so- 
lution. In two of the three examples problems consid- 
ered. however, using the initial vectors in an overall 
finite difference method or assuming fixed modes in 
a semianalytical method resulted in very poor modal 
convergence for stress sensitivities. Two methods 
were devised to improve this poor performance. 

The first method retains the derivatives of the 
basis vectors in the sensitivity equations but approx- 
imates these derivatives rather than using a very 
costly exact computation. One well-known method 
for approximating eigenvector derivatives, the modal 
method, was found to be completely ineffective be- 
cause it adds no new information to the existing 
modal basis. Another technique, the modified modal 
method, adds a pseudostatic contribution to the 
eigenvectors in approximating the eigenvector deriva- 
tives. This technique, along with the semianalytical 
method, was found to be very effective in improving 
the poor accuracy of the stress sensitivities. 

A second method for improving the accuracy of 
the stress sensitivities as a function of the number of 
modes is to use a mode acceleration version of the 
semianalytical method. The key to the mode accel- 
eration method in the transient analysis is that it 
supplements the modal basis with a static contribu- 
tion calculated from the complete model. The key to 
the mode acceleration implementation of the semi- 
analytical sensitivity method is that it supplements 
the modal basis with pseudostatic sensitivity terms 
calculated from the complete model. This technique 
produced the same dramatic improvement in the ac- 
curacy of stress sensitivities as the semianalytical 
modified modal method. 

As mentioned, computational cost was an over- 
riding concern in considering the sensitivity analy- 
sis methods. To estimate this cost, expressions for 
the number of floating point operations in each of 
the methods were derived. Although this approach 
does not include important effects such as overhead 
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operations or disk input/output that would be 
present in a practical implementation of these meth- 
ods, it does provide a mechanism for an approxi- 
mate coarse ranking of the methods by computa- 
tional cost. The overall forward difference method 
with fixed basis vectors was found to be the cheap- 
est method for all cases considered. This technique, 
however, suffers from the accuracy problems previ- 
ously mentioned. One approach to alleviating these 
accuracy problems is to recalculate the modes for 
the perturbed model (updated modes) in the overall 
forward difference method. This forward difference 
method with updated modes was found to be very 
costly for large models, however. The fixed-mode 
semianalytical method is only slightly more costly 
than the overall forward difference method with fixed 
modes but suffers from the same accuracy problem 
as the fixed-mode overall forward difference method. 
Two techniques with reasonable costs that alleviate 
the accuracy problem are the mode acceleration im- 
plementation of the semianalytical method and the 
semianalytical method with approximate mode shape 
derivatives. Of these two methods, the semianalyti- 
cal mode acceleration method is slightly cheaper. 

Given the high accuracy of the semianalytical 
mode acceleration method for a relatively small num- 
ber of modes and its reasonable computational cost, 
this appears to be the method of choice. In the three 
examples considered herein, this method consistently 
performed as well as the much more costly, updated- 
mode overall finite difference methods. Furthermore, 
the insensitivity of this and the other semianalyti- 
cal methods to finite difference step size makes this 
semianalytical mode acceleration method especially 
attractive. 
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Appendix 

Computer Implementation 


The methods for calculating sensitivities and the example problems have been implemented 
with the general purpose finite clement code, EAL (ref. 23). EAL includes general language 
constructs for controlling execution flow as well as general and specific utilities for manipulating 
data stored as named entities in a data base. It also allows procedures (calle runstreams ) 
to be defined and then explicitly executed. Most of the implementation was done with EAL 
runstreams. However, some parts of the implementation could not be conveniently done va 
runstreams and were coded as Fortran additions to EAL. The Fortran additions are descn >e 
in the next section. The runstreams for the algorithms and example problems are included an 

described also. 


Additions to EAL 

The transient response module in EAL version 312 solves the uncoupled form of equa- 
tions (2.5) with the matrix series expansion method. A modification was made to allow equa- 
tions (2.5) to be fully coupled. In the scmianalytical method, the right-hand-side pseudo loading 
of equations (4.6) can be easily formed with EAL. However, a slight modification to the tran- 
sient response module was required to permit solution of equations (4.6) with this general form 
of loading. In addition, a special purpose module was added to EAL to perform the task of 
identifying the critical points on each response function. 


Runstream for Stepped Beam Example Problem 

The runstream for the stepped beam is included to illustrate how the sensitivity calculation 
runstreams are used. At the beginning of the runstream, the data set XFLG ADS 'ideates whic i 
subset of the possible design variables will be considered in the sensitivity analysis. The data 
sets X ADS and XNAME ADS contain the initial values and the register names of all the design 
variables, respectively. Various parameters controlling the analysis and sensitivity calculations 
are defined in runstream data sets TR PARAMETERS. DXDV PARAMETERS, and BACK METHOD. The 
runstream data set MODEL defines the model in terms of the design variables in X ADS It is cal e 
before the initial dynamic analysis and at least once for each design variable considered in the 
sensitivity analysis. The runstream data set DYNAM SOLN is called once to perform the dynamic 
analysis of the initial model. The runstream data set PLOT RESP illustrates the interface to a 
useful utility runstream TR PLOT for automatically generating plots of response quantities as a 
function of time. TR PLOT is called once for each class (e.g., accelerations) of response quantity 
to be plotted. The actual sensitivity analysis is performed by calling the runstream TR DXDV n 
where the n is associated with the particular sensitivity calculation method. 


+CM=120000 

$ 

$ RUNSTREAM FOR STEPPED BEAM EXAMPLE 

$ 

*XQT EXTE 

! SYST = SSP(4 ,5) $ GET SYSTEM TYPE 
*XQT U1 



*INF=7 
*CLIB=29 
'* (ALL) ALL 
*XQT AUS 

TABLE(NI=1 , NJ=9 , TYPE=0) : XFLG ADS 
J=1 : 1 
J=5 : 1 
J=6 : 1 
J=9 : 1 
*XQT U1 
*TI(X ADS) 

23.5 

22.0 

20.0 

18.0 

16.5 

$ 

40.0 

80.0 

120.0 

160.0 

♦TKXNAME ADS) 

HI : H2 : H3 : H4 : H5 
XL1 : XL2 : XL3 : XL4 
* (TR , PARAMETERS) 

QLIB=1 
MNAME=CEM 
NM0DES=5 
DT-1 . OE-5 


T2 = .3 
DRF0RMAT=DIAG 
METHOD-MODES 
DXDV=0 
EIGEN=1 
PRINT=1 
VLIB-1 
NCRIT=5 
C0NV-1.E-10 
BLKSIZE=2000 
NTERMS=50 
* ( DXDV , PARAMETERS ) 
FDCH=1 . OE-5 


FDMCH=1 .OE-6 


DXMD=FIXED 


*TI (BACK METHOD) 

2 $ DISP 
1 $ VELOCITIES 
1 $ ACCELERATIONS 

1 $ REACTIONS 

2 $ STRESSES 

* (MODEL) E ND 

! LEN=200 . 

! NEPS = 3 
! NEL » NEPS *5 
! NNODE * NEL + 1 


*XQT AUS 

TABLE(NI=1,NJ=5) : XXI 

1=1 : J=1 , 5 : 0. "XL1 " "XL2" "XL3" "XL4" 

TABLE(NI=l,NJ-5) : XX2 

1=1 : J=l,5 : "XL1 " ”XL2" "XL3" "XL4" "LEN" 
D1 = SUM(XX2 -1.0 XXI) 

! RNEP = 1.0/NEPS 


DELX = UNION ("RNEP" Dl) 
*XQT TAB 

START "NNODE" 1345 
JLOC 
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$10. 0. 0. 200. 0. 0. "NNODE" 

! J = 1 : ! X = 0.0 
! II = 1 : ! N1 = 5 
♦LABEL 5 

! DELX = DS,1 ,"I1" ,1(1 DELX AUS 1 1) 

! 12 = 1 : ! N2 = NEPS 
♦LABEL 8 
’’J” "X" 0.0 0.0 
! J = J + 1 : ! X = X + DELX 
♦ JGZ , -1 (N2 , 8) 

! II = II + 1 
* JGZ, -1 (N1 , 5) 

"NNOD" "LEN" 0. 0. 

CON 1 

ZERO 1,2,6 : 1 
MATC 

1 30. +6 .3 .3 
BA 

RECT 1 1.20 "HI" 

RECT 2 1.10 "H2" 

RECT 3 1.00 "H3" 

RECT 4 .90 "H4" 

RECT 5 .85 "H5" 

RECT 6 1.0 20.0 
MREF 

1 1 2 1 1.0 
♦XQT ELD 
E21 

! N = 5 
! I = 1 
! N1 = 1 
♦LABEL 20 
NSECT = "I" 

! N2 = N1 + 1 
"Nl" "N2" 1 "NEPS" 

! Nl = Nl + NEPS 
! 1 = 1 + 1 
♦JGZ, -1 (N ,20) 

♦XQT E 
RESET G=386 . 

♦XQT EKS 
♦XQT TAN 
♦XQT K 
♦ XQT M 
RESET G=386 . 

♦XQT AUS 
R = RIGID(l) 

DEFINE R6 = R AUS 1 1 6,6 
APPL FORC = PROD (-1.0 CEM R6) 

♦END 

♦ (DYNAM , SOLN) END 
♦DCALL(TR, VECTORS) 

♦XQT U1 
♦TICSEL DISP) 

"NNODE" 2 
♦TICSEL VELO) 

"NNODE" 2 
♦TICSEL ACCE) 

"NNODE" 2 
♦TICSEL STRE) 

E21 1 1 SZ1 1 0 
♦XQT AUS 

$ DEFINE SOME MODAL DAMPING 
TABLE (NI=1 ,NJ="NMODES") : DRAT 
1=1 : J=1 , "NMODES" : .005 
♦DCALL (SQUARE LOAD) RTIME=.18 RANG=10.0 


♦DCALL ( TR , MA IN ) 


END 


♦(SQUARE, LOAD) END 
♦XQT AUS 

! RT2 = RTIME/2 . 0 
! EPS = 0.0 
! RT2M = RT2 - EPS 
* RT2P = RT2 + EPS 
! RTM = RTIME - EPS 
! D2R = 3.1415926/180. 

! RRAD = RANG*D2R 
! AMP = 4 . 0*RRAD/RTIME/RTIME 
! MAMP = -AMP 
TABLE(NJ=6) : TIME 

J=l,6 : 0. ,, RT2M ,, "RT2P" "RTM" "RTIME" 10000.0 
TABLE(NJ=6) : CA 

J=1 ,6 : "AMP" "AMP" "MAMP" "MAMP" 0.0 0.0 
♦END 

♦(PLOT RESP) END 
♦XQT DCU 

CHANGE 1 A AUS MASK MASK HIST CA 1 1 
♦XQT AUS 
ALPHA : FTITLE 

1 ’ HISTORY OF FORCE MULTIPLIER G(T) 

ALPHA : DTITLE 

1* TIP DISPLACEMENT HISTORY FOR CANTILEVER BEAM 
ALPHA : TVTITLE 

1 ’ TIP VELOCITY HISTORY FOR CANTILEVER BEAM 
ALPHA : TATITLE 

1’ TIP ACCELERATION HISTORY FOR CANTILEVER BEAM 
ALPHA : STITLE 

1’ BENDING STRESS AT THE ROOT FOR THE CANTILEVER BEAM 
! TLIB = 15 
♦XQT U1 

♦(TRPLOT OPTIONS) 

! YNAME = ’CA 
! TITLE = ’FTITLE 
! ID = 1 
♦DCALL(TR.PLOT) 

♦XQT U1 

♦(TRPLOT OPTIONS) 

! YNAME = ’DISP 
! IDJK = ’DISP 
! TITLE = ’DTITLE 
! ID = "NMODES" 

♦DCALL(TR,PLOT) 

♦xqT U1 

♦(TRPLOT OPTIONS) 

! YNAME = ’VELO 
! IDJK = ’VELO 
! TITLE = ’TVTITLE 
! ID = "NMODES" 

♦DCALL(TR.PLOT) 

♦XQT U1 

♦(TRPLOT OPTIONS) 

! YNAME = ’ACCE 
! IDJK = ’ACCE 
! TITLE = ’TATITLE 
! ID = "NMODES" 

♦DCALL ( TR , PLOT ) 

♦xqT U1 

♦(TRPLOT OPTIONS) 

! YNAME = ’STRE 
! IDQ = ’STRE 
! TITLE * ’STITLE 
! ID = "NMODES" 
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♦DCALL (TR, PLOT) 

♦END 

♦RGI 

♦DCALL (TR , PARAMETERS ) 

♦DCALL(SENS,DVUP) 

♦DCALL (MODEL) 

♦DCALL (DYNAM , SOLN) 

$^DCALL(PLOT,RESP) 

♦IFC’DXDV" NE 0): +DCALL (TR , DXDV , "DXDV" ) 

♦ ALL 

♦ IFC’SYST" EQ CDC ): ♦PRT(ALL) 

♦ IFC'SYST" Eq CNVX) : ♦PRT(ALL) 

♦DCALL(ALL) 

♦xqT EXIT 

Runstreams for Sensitivity Methods 

Runstream TR MAIN 

TR MAIN is the main runstream for performing the transient response analysis and is based 
on a similar runstream produced by EISI. It is used only for the transient analysis of the initial 
model and not for the sensitivity calculations. 


$ (TR MAIN) - MAIN DRIVER FOR TRANSIENT RESPONSE ANALYSIS 

$ 

♦ XQT U1 

♦ REGISTER STORE (29 TR REGISTERS 1 1) 

♦ REGISTER RETR (29 TR REGISTERS 1 1) 

♦ RGI 

$ DEFAULT REGISTERS: „„„„„ 

QLIB =2 $ SOURCE FOR EXCITATION, DESTINATION LIB FOR RESPONSE 

VLIB =1 $ SOURCE LIB FOR VIBR MODE AND VIBR EVAL DATASETS 

VSET =1 $ USE "VLIB" VIBE MODE "VSET" "VCON" FOR THE RITZ 

VCON * 1 $ FUNCTIONS 

KNAME = K $ STIFFNESS MATRIX 

MNAME = DEM $ MASS MATRIX 

DAMP = DAMP $ NAME OF SPAR FORMAT DAMPING MATRIX 

FSET * 1 $ EXCITATION SET NUMBER 

NAME = AUS $ 2ND WORD OF TIME "NAME" AND CA "NAME" 

DT =0. $ SET TIME INCREMENT 

T2 =0. $ FINAL INTEGRATION TIME 

DRFORM= DIAG $ FORMAT FDR THE REDUCED MATRICES (DIAG .FULL , RITZ) 
DRMETH= 0 $ TIME INTEGRATION METHOD (0=MSE) 

NTER =50 $ SET NUMBER OF TERMS IN MATRIX SERIES EXPANSION 

NMODES* 0 $ NUMBER OF MODES USED IN DYNAMIC ANALYSIS (DEFAULT* ALL) 

BLKSIZ= 6000 $ BLOCK SIZE FOR OUTPUT DATASETS 

EIGEN =0 $ EIGENVALUE ANALYSIS OF DAMPED SYSTEM 

PRINT =0 $ PRINT FLAG FOR DTEX 

OPT=0 , PROC=MAIN , NERR=0 

♦ DCALL, OPT (TR PARAMETERS) 

» ZERO = SSP(0 , 10) 

♦IFC'NMODES" Eq 0): ! NMODE=TOC , NBLOCK( "VLIB" VIBR MODE "VSET" "VCON") 

$ 

$ COMPUTE DATASETS REQUIRED FOR DR/DTEX , /TR1 , AND /BACK : 

$ 

♦ CALL (TR PREP) 

$ 

$ COMPUTE THE MODAL RESPONSE: 

$ 

♦ XQT DRX 

♦ IFC’DT" GT l.E-20): ♦GOTO 20 

DTEX(INLI="QLIB" ,N2="NAME" , OUTL="QLIB" , EIGEN="EIGEN" ,> 

PRINT="PRINT" ) 


♦GOTO 30 



♦LABEL 20 


♦LABEL 30 


DTEX(INLI="QLIB" ,N2="NAME H »0UTL="QLIB" ,DT="DT" 
NTER= 11 NTER" , EIGEN="EIGEN" , PRINT= "PRINT " ) 


> 


TR1 ( INLIB="QLIB" ,N2="NAME M f CASE="FSET" , ALIB="QLIB"> 
QXLIB="QLIB" , QXlLIB="qLIB" ,QX2LIB="QLIB" ,> 

T2="T2" , LB="BLKSIZ") 


$ 

$ BA%£ TRANSFORMATION: 

♦XQT MJS 
! NBCK = 5 

TABLE(NI="NBCK\ NJ=1, TYPE=4) : "QLIB" BACK LIST 
J=1 : DISP VELO ACCE REAC STRE 


$ 


EXIT 

! I = 1 : ! N = NBCK 
♦LABEL 50 

• BKMETH = DS , 1 , " I" , 1 ("qLIB" BACK METH 1 1 ) 

! NM = DS/'I'M.lCqLIB" BACK LIST 1 1) 

! IERR * TOC, IERRO'QLIB" SEL "NM" MASK MASK) 
♦IFO'IERR" EQ 0): ♦CALL (TR "NM" "BKMETH") 

! I = I + 1 

*JGZ,-1(N,50) 

$ 

$ EXIT: 


♦ XQT U1 

♦ REGISTER RETRIEVE (29 TR REGISTERS 1 1) 
♦END 


Runstream TR PREP 


This Runstream TR PREP is used to define the reduced equations for the transient analysis 
and also to prepare data for the back transformation phase. It is based on an EISI runstream 
of the same name. It is used only for the transient analysis of the initial model and not in the 
sensitivity calculations. 


$ 

$ (TR PREP) - PREPARATION OF REQUIRED DATASETS 
$ 

♦ XQTC U1 

♦ RGI 

PRDC=PREP , NERR=0 , M0TI=M0TI, F0RC=F0RC 

$ 

$ CHECK FDR THE REQUIRED DATASETS AND DETERMINE THE TYPE OF EXCITATION: 
$ ! TYPE=0 FOR APPLIED FORCE. !TYPE=1 FOR APPLIED DISPLACEMENT 
$ 


! TYPE=2 

! IERR=T0C , IERR( "QLIB" APPL FORC "FSET" MASK) 

♦ IF ( "IERR" NE 0): ♦GO TO 100 

! TYPE-0 : !FNAME= J F0RC $$ APPLIED FORCE EXCITATION 

♦ LABEL 100 


!IERR=T0C, IERRO'QLIB" APPL MOTI "FSET" MASK) 

♦ IF ("IERR" NE 0): ♦GO TO 200 

♦ IF ( "TYPE" EQ 0): ! NERR= 1 $$ FORCE & DISP SPECIFIED 

! TYPE= 1 : !FNAME= , MOTI $$ APPLIED DISPLACEMENT EXCITATION 

♦ LABEL 200 


♦ IF ("TYPE" EQ 2) : ! NERR=2 j 

♦ IF ( "NERR" NE 0) : ♦CALL (TR ERROR) 

! IERR=T0C , IERR ( " VLIB 11 VIE 

♦ IFO'IERR" NE 0) : ♦CALL (TR ERROR) 

! IERR=T0C , IERR ( "VLIB" VIE 

♦ IF ("IERR" NE 0) : ♦CALL (TR ERROR) 

! IERR=TDC, IERRO'QLIB" TIN 

♦ IF ("IERR" NE 0) : ♦CALL (TR ERROR) 


NO EXCITATION SPECIFIED 

$ NERR=2 

MODE ' 

’VSET" "VCON") : 

! NERR=3 

$ NERR-3 

EVAL ' 

•VSET" "VCON"): 

! NERR=4 

$ NERR =4 

"NAME 1 

' "FSET" MASK): 

! NERR=5 

$NERR=5 
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! IERR=TOC , IERR ( "QLIB" CA 
IF ("IERR" NE 0) : *CALL (TR ERROR) 


NAME" "FSET" MASK): ! NERR=6 


$NERR=6 


$ 

$ COMPUTE XTMX , XTKX t XTDX, & XTF FOR DR/DTEX & TR1 : 

$ 

! IERR=TOC , IERR ( 1 INV "KNAME" "VCON" MASK) 

♦ IF ("IERR" EQ 0) : *G0 TO 250 

♦ XQT DRSI 

RESET K="KNAME" , CON="VCDN" 

♦ LABEL 250 

$ 

*XQT AUS 


0UTLIB= 

="QLIB" : 

INLIB=' 

"QLIB" 



DEFINE 

X="VLIB" 

VIBR 

MODE 

"VSET" 

"VCON 1 

DEFINE 

E=”VLIB” 

VIBR 

EVAL 

"VSET" 

"VCON 

DEFINE 

F="QLIB" 

APPL 

"FNAME" 

"FSET" 



DEFINE K= 1 "KNAME” 

DEFINE M= 1 "MNAME” 

♦ IF( "TYPE" EQ 0): *G0 TO 300 

KS=PROD("KNAME" -1. F) : DEFINE F=KS 

♦ LABEL 300 

XTF "NAME” "FSET"= XTY(X.F) 

$ 

! IDMD = TOC , IERR( "VLIB” DRAT MASK MASK MASK) 

♦ IFC’IDMD" NE 0): ♦GOTO 400 

DEFINE D = "VLIB” DRAT 
OMEG = SQRT(E) 

DMPD = PROD (2.0 D OMEG) 

♦LABEL 400 

! IERR=TOC , IERRO'QLIB" XTMX "NAME” MASK MASK) 

♦ IF ("IERR" EQ 0): ♦GO TO 500 

$ 

*CALL(TR,REDM) 

♦ LABEL 500 

$ 

$ COMPUTE ("QLIB” STAT DISP "FSET" "VCON"): 

$ 

! I ERR=TOC, IERRC’QLIB" STAT DISP "FSET" "VCON") 

♦ IFO'IERR" EQ 0): *G0 TO 700 

♦ XQT SSOL 

RESET SET=”FSET" ,CON="VCQN" ,QLIB="QLIB” 

♦ LABEL 700 

$ 

♦LABEL 1000 
♦END 

Runstream TR DISP 

The two TR DISP runstreams do the back transformation for displacements. TR DISP 1 
performs the back transformation with the mode displacement method and TR DISP 2 uses 
the mode acceleration method. This naming convention is used for the other runstreams 
that perform the back transformation operation for other response quantities. The TR DISP 
runstreams and companion runstreams for velocities, accelerations, and stresses perform the 
back transformation at all time steps. Accordingly they are used only in the dynamic analysis 
of the initial model and not in the sensitivity analysis. 

$ 

$ (TR DISP 1) - BACK TRANSFORMATION FOR DISPLACEMENTS 
$ MODE DISPLACEMENT METHOD 

$ 

♦XQTC AUS 

OUTLIB="QLIB” : INLIB="QLIB" 

DEFINE IDJK = "QLIB" SEL DISP 

DEFINE X = "VLIB" VIBR MODE "VSET" "VCON" 


TMAT VMOD=SVTRAN(IDJK ,X) 

♦XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1.0 "QLIB" TMAT VMOD : Y = "QLIB" QX 

Z~ "QLIB" HIST DISP 

EXT = "QLIB" EXT DISP "FSET" 

♦END 

$ 

$ (TR DISP 2) - BACK TRANSFORMATION FOR JOINT DISPLACEMENTS 
$ MODE ACCELERATION METHOD 

$ 

$ 

$ TRANSIENT RESPONSE: BACK TRANSFORMATION FOR JOINT DISPLACEMENTS 
$ APPLICABLE FOR APPLIED FORCE OR DISPLACEMENT EXCITATION 
$ 

$ REGISTERS: QLIB, NAME, FSET, VCON 
$ 

♦ XQT AUS 

! ID = TOC, IERR(" QLIB" XTDX MASK MASK MASK) 

OUTLIB="QLIB" : INLIB="QLIB" 

DEFINE E = "VLIB" VIBR EVAL 
ROMG = RECIP(E) 

DEFINE IDJK = "QLIB" SEL DISP 

DEFINE XS = "QLIB" STAT DISP "FSET" "VCON" 

DEFINE X = "VLIB" VIBR MODE 1 "VCON" 

DEFINE DACC = TMAT DACC 

XOME = CBD (X , ROMG) 

TMAT DACC = SVTRAN(IDJK ,XOME) 

TMAT DS = SVTRAN(IDJK.XS) 

♦JNZ(ID ,30) 

$ DAMPING TERM 

! NJ = TOC,NJ ( "QLIB" XTDX MASK MASK MASK) 

! NINJ = TOC,NINJ ("QLIB" XTDX MASK MASK MASK) 

*IF( "NJ" NE "NINJ"): *GOTO 10 
$ MODAL DAMPING 

XOMD = CBD (XOME, XTDX) 

TMAT DVEL = SVTRAN ( IDJK , XOMD) 

♦GOTO 30 
♦LABEL 10 
$ GENERAL DAMPING 

TMAT DVEL = RPROD (DACC , XTDX) 

♦LABEL 30 
♦ XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1. "QLIB" TMAT DS : Y = "QLIB" A "NAME" "FSET" 
♦IF( " ID" EQ 0): T*-l. "QLIB" TMAT DVEL : Y="QLIB" QX1 "NAME" "FSET" 


T = -1. 

"QLIB" 

TMAT 

DACC 

: Y = "qLIB" QX2 "NAME" "FSET" 

Z 

"QLIB" 

HIST 

DISP 

"FSET" 

EXT = 

"QLIB" 

EXT 

DISP 

"FSET" 


♦END 

Runstream TR VELO 


$ 

$ (TR VELO 1) - BACK TRANSFORMATION FOR VELOCITIES 
$ MODE DISPLACEMENT METHOD 

$ 

♦XQTC AUS 

OUTLIB="QLIB": INLIB="QLIB" 

DEFINE IDJK = "qLIB" SEL VELO 

DEFINE X * "VLIB" VIBR MODE "VSET" "VCON" 
TMAT VVEL=SVTRAN(IDJK,X) 

♦XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1.0 "QLIB" TMAT VVEL : Y = "QLIB" QX1 

66 


♦END 


Z= "QLIB" HIST VELO 

EXT = "QLIB" EXT VELO "FSET" 


Runstream TR ACCE 


$ 

$ (TR ACCE 1) - BACK TRANSFORMATION FOR ACCELERATIONS 
$ MODE DISPLACEMENT METHOD 

$ 

♦ XQTC AUS 

OUTLIB="QLIB" : INLIB="QLIB" 

DEFINE IDJK = "QLIB" SEL ACCE 

DEFINE X = "VLIB" VIBR MODE "VSET" "VCON" 
TMAT VACC-SVTRAN (IDJK , X) 

♦XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1.0 '’QLIB" TMAT VACC : Y = ’'QLIB 11 QX2 

Z= "QLIB" HIST ACCE 

EXT = "QLIB" EXT ACCE "FSET" 

♦END 

Runstream TR STRESS 


$ 

$ (TR STRESS l) 

$ 

$ 

$ MODE DISPLACEMENT STRESS BACK TRANSFORMATION 

$ 

♦XQT ES 
RESET OPER=T 
IDQ= "QLIB" SEL STRESS 
U = "VLIB" VIBR MODE 1 "VCON" 1 , "NMDDE" 

T = "QLIB" TMAT VSTRE "FSET" 

♦XQT DRX 

BACK (LRZ= s "BLKSIZE") 

T = +1. "QLIB" TMAT VSTRE "FSET" : Y = "QLIB" qX "NAME" "FSET" 
Z = "QLIB" HIST STRESS "FSET" 

EXT *= "QLIB" EXT STRESS "FSET" 

♦END 

$ 

$ (TR STRESS 2) 

$ 

$ 

$ TRANSIENT RESPONSE: BACK TRANSFORM FDR ELEMENT STRESSES 
$ APPLICABLE FOR APPLIED FORCE OR DISPLACEMENT EXCITATION 
$ 

$ REGISTERS: QLIB, NAME, FSET, VCON 

$ 

♦ XQT AUS 

! ID = TOC , IERRC'QLIB" XTDX MASK MASK MASK) 

OUTLIB="QLIB" : INLIB="QLIB" 

DEFINE E = "VLIB" VIBR EVAL 
ROMG = RECIP(E) 

DEFINE IDJK = "QLIB" SEL DISP 

DEFINE XS = "QLIB" STAT DISP "FSET" "VCON" 

DEFINE X = "VLIB" VIBR MODE 1 "VCON" 

XOME = CBD(X ,ROMG) 

*JNZ(ID ,30) 

$ DAMPING TERM 

! NJ = TOC ,NJ ("QLIB" XTDX MASK MASK MASK) 

! NINJ = TOC , NINJ ("QLIB" XTDX MASK MASK MASK) 

♦ IFC'NJ" NE "NINJ"): ♦GOTO 10 
$ MODAL DAMPING 
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XOMD = CBD ( XOME , XTDX ) 


♦GOTO 30 
♦LABEL 10 
$ GENERAL DAMPING 

XOMD=CBR (XOME, XTDX) 

♦LABEL 30 
♦XQT ES 

RESET OPER=T 

IDQ = "QLIB" SEL STRESS 

u= M q 

♦IF("ID" 

♦ XQT DRX 


♦ IFC’ID" 


U= "qLIB" STAT DISP "FSET" 
EQ 0): U= "QLIB" XOMD : 

U= "QLIB" XOME : 

: 

"VCON" 

: T= » 
T= 
T= 

QLIB" TMAT SF 
"QLIB" TMAT SD 
"QLIB" TMAT SP 

BACK(LRZ="BLKSIZE") 

T = +1. "qLIB" TMAT SF 

: Y = 

"QLIB" 

A "NAME" "FSET" 

EQ 0): T— 1. "QLIB" TMAT SD : 

Y = 

"QLIB" 

QX1 "NAME" "FSET" 

T * -1 . "QLIB" TMAT SP 

: Y = 

"QLIB" 

qX2 "NAME" "FSET" 


♦END 


Z * "qLIB" HIST STRE "FSET" 

EXT = "qLIB" EXT STRE "FSET" 


Runstream TR ERROR 


$ 

$ (TR ERROR) 

$ 

$ 

$ THIS PROCEDURE PRINTS FATAL ERROR MESSAGES FOR THE TR PROCS 

$ 

♦ XqT U3 

RP2 : NUMBER OF F0RMATS*10 

FORM l’(33Hl*** TR FATAL ERROR: PROC, NERR= ,A4,1H,,I4) 

PRINT(l) "PROC" "NERR" 

♦ GO TO "PROC" 

$ 

♦ LABEL PREP 

$ 

FORM l'ClOX^/H .BOTH (APPL FORC FSET ) AND (APPL MOTI FSET )/ 

10X , 46HARE SPECIFIED IN qLIB , ONLY ONE IS PERMITTED) 

FORM 2’ ( 10X,48HNEITHER (APPL FORC FSET ) NOR (APPL MOTI FSET )/ 

* 10X , 20HIS PRESENT IN QLIB ) 

FORM 3 ’ (10X ,47HVIBR MODE VSET VCON IS NOT PRESENT IN VLIB ) 

FORM 4’ (10X .47HVIBR EVAL VSET VCON IS NOT PRESENT IN VLIB ) 

FORM 5 ’ ( 10X , 43HTIME NAME FSET IS NOT PRESENT IN QLIB ) 

FORM 6’ (10X.43HCA NAME FSET IS NOT PRESENT IN QLIB ) 

PRINT ( "NERR") 

$ 

♦ GO TO FINIS 

$ 

♦ LABEL FINIS 

♦ XqT U1 

♦ SHOW 

♦ XQT DCU 

TOC 1: TOC "QLIB": TOC "VLIB" 

♦ XQT EXIT 
♦END 

Runstream TR RITZ 

Runstream TR RITZ calculates RWL vectors following equations (2.25) through (2.29). Then 
the reduced system is optionally uncoupled by solving the reduced-order eigenproblem. This 
runstream is substantially based on one written at EISI. 

$ 

$ (TR RITZ) 

$ 
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* XQT U1 
$ 

$$$ REGISTERS: MAXHZ , AFLIB, VLIB, MNAME , NMD, SCALE 

$ 

* REGISTER STORE (29 REGISTER HOUSE 1 1) 

* REGISTER RETR (29 REGISTER HOUSE 1 1) 

! IERR=TOC IERR(1 INV K MASK MASK) 

*ONLINE=0 

* IFC'IERR" EQ 0): *G0 TO 109 

* XQT DRSI 

* LABEL 109 

* XQT SSOL 

RESET QLIB-’ 1 AFLIB" 

* XQT AUS 

0UTLIB=10: INLIB=10 
DEFINE M=1 "MNAME" 

$ 

$ SCALE THE FIRST VECTOR 

$ 

DEFINE X-"AFLIB" STAT DISP 1 MASK 1 
MX=PROD(M,X) 

XTMX=XTY (X , MX) 

RECI=RECI (XTMX) 

SCAL=SQRT(RECI) 

11 RITZ VECT =CBD ( X , SCAL ) 

DEFI RITZ=11 RITZ VECT 

12 MX=PROD(M ,RITZ) 

! NSET=TDC NBLOCKS ( "AFLIB" STAT DISP 1 MASK) 

!N=NSET-1: !N1*1 : ! N2=0 

* IF ("N" Eq 0): *G0 TO 104 

$ 

$ M-ORTHONORMALIZE VECTORS 2 THROUGH NSET 

$ 

* LABEL 105 

tNl-Nl+1 : ! N2-N2+1 

DEFI U="AFLIB" STAT DISP 1 MASK "Nl" 

DEFI MX*12 MX MASK MASK MASK 1 "N2" 

XTMU=XTY (MX ,U) 

DEFI X-ll RITZ VECT MASK MASK 1 "N2" 

A=CBR(X , XTMU) 

U1=SUM(U, -1 . A) 

MU=PR0D(M,U1) 

UTMU=XTYD(U1 ,MU) 

REC I =REC I ( UTMU ) 

SCAL=SQRT (RECI ) 

VECT=CBD (U1 , SCAL) 

11 RITZ VECT=UNION ,U(VECT) 

TEMP^PROD ( M , VECT ) 

12 MX =UNION ,U(TEMP) 

* JGZ -1 (N 105) 

* LABEL 104 

* XQT DCU 

ERASE 10 

! NDO-NMD-NSET : »SET*NSET+1: !SET1=NSET 

* IF ("NDO" LE 0): ♦GO TO 1002 

* XQT AUS 

TABLE(NJ=1) : 13 SCALE 
J=1 : "SCALE" 

$ 

$ GENERATE REMAINING VECTORS ORTHONORMAL TO STATIC SOLUTION RITZ VECTORS 
$ 

* LABEL 1000 

* XQT AUS 

0UTLIB=10 : INLIB=10 
DEFINE M=1 "MNAME" 

DEFINE X=ll RITZ VECT MASK MASK "SET1" 
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TEMP-PROD (M,X) 

NORM-NORM (TEMP) 

DEFI SCALE-13 SCALE 
APPL FORC-CBD (NORM , SCALE) 

* IF ( "SET1" EQ "NSET"): *GO TO 106 
12 MX-UNION , U(TEMP) 

* LABEL 106 

* XQT SSOL 

RESET QLIB-10 

* xqT AUS 

0UTLIB-10 : INLIB-10 
DEFI U-STAT DISP 

DEFI MX-12 MX MASK MASK MASK 1 "SET1" 
XTMU-XTY (MX ,U) 

DEFI X=ll RITZ VECT MASK MASK 1 "SET1" 
A=CBR(X , XTMU) 

U1=SUM(U, -1 . A) 

DEFI M-l "MNAME" 

MU=PROD(M, Ul) 

UTMU-XTYD (Ul , MU) 

RECI-RECI (UTMU) 

SCAL=SQRT(RECI ) 

VECT=CBD (U 1 , SC AL ) 

11 RITZ VECT=UNION,U(VECT) 

! SET1-SET : !SET=SET+1 

* XQT DCU 

ERASE 10 

* JGZ.-KNDO, 1000) 

* LABEL 1002 

♦IF("DRFORMAT" EQ DIAG) : ♦ GOTO 10020 

* XQT AUS 

DEFI X=ll RITZ VECT 

"VLIB" VIBR MODE 1 1=UNI0N(X) 

TABLE (NI-1 , NJ-"NMD" ) : VIBR EVAL 
* RETURN 
♦LABEL 10020 

* XQT AUS 

0UTLIB-10 : INLIB-10 

DEFINE K=1 K: DEFI M=1 "MNAME" 

DEFINE X=ll RITZ VECT 
I JC0DE=10000 
! NMODE-NMD 

KX=PRDD(K ,X) : SYN K 10000 "NMODE" = XTYS(X.KX) 
MX=PROD(M ,X) : SYN M 10000 "NMODE" = XTYS(X,MX) 
! ZERO=NMODE- 1 

* JZ (ZERO, 1003) 

♦XQT DCU 

TOC 10 

* XQT STRP 

RESET SOURCE- 10 , DEST-10, FRQ2="MAXHZ" 

* JGZ (ZERO, 1004) 

* LABEL 1003 

* XqT AUS 

OUTLIB-IO: INLIB-10 
• K-DS 2 1 1(10 SYN K MASK MASK) 

• M=DS 2 1 1(10 SYN M MASK MASK) 

! EVAL-K/M 

TABLE(NI=1 ,NJ=1) : SYS EVEC : J=l: 1.0 
TABLE(NI=1 ,NJ=1) : SYS EVAL: J=1 : "EVAL" 

* LABEL 1004 

* XQT AUS 

OUTLIB-IO: INLIB-10 
DEFINE E-SYS EVEC 
DEFI X-ll RITZ VECT 
X ORTH 1 1=CBR(X,E) 

DEFINE X-X ORTH 1 1 
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"VLIB" VIBR MODE 1 l=UNION(X) 

DEFINE E=SYS EVAL 

"VLIB" VIBR EVAL 1 l=UNION(E) 

"VLIB" VIBR HZ 1 1=SQRT( .0253303 E) 
♦0NLINE=1 

* XQT DCU 

PRINT "VLIB" VIBR HZ 1 1 

* XQT U1 

* REGISTER RETR (29 REGISTER HOUSE 1 1) 

♦END 


Runstream TR REDM 

TR REDM is a utility runstream for generating the reduced equations given a set of basis 
vectors. Depending on the input register DRFORMAT the equations can be coupled or uncoupled. 
If the equations are uncoupled, it is assumed that is the identity matrix and <I> K4> 

is a diagonal matrix with the eigenvalues along the diagonal. 

$ 

$ (TR REDM) - FORM REDUCED K AND M MATRICES FOR TRANSIENT RESP . 

$ 

$ 

$ REGISTERS: 

$ DRFO = 'FULL, ’RITZ, OR 'DIAG 

$ NMODE * NUMBER OF MODES 

$ MNAME = MASS MATRIX NAME 

$ VLIB = LIBRARY FDR VIBRATIONAL MODES AND FREQS 

$ QLIB = DESTINATION LIBRARY FOR MATRICES 

$ 

♦ XQTC AUS 
OUTLIB="QLIB" 

DEFINE X = "VLIB" VIBR MODE 1 1 1, "NMODES" 

DEFINE E = "VLIB" VIBR EVAL 1 1 

DEFINE DAMP = 1 DAMP SPAR 
DEFINE DMPD = 1 DMPD 

! IDSP = TOC, IERR(1 DAMP SPAR MASK MASK) 

! IDMD = TDC , IERR( 1 DMPD MASK MASK MASK) 

! DRFO 

♦IF("DRFO" NE FULL): ♦GOTO 100 
$ 

$ FULL MATRICES, X IS A SET OF EIGENVECTORS 

$ 

! N = NMODES 
! I = 1 

TABLE(NI="NMODES" ,NJ="NMODES" ) : XTMX 
♦LABEL 10 

I="I" : J="I" : 1.0 
! 1 = 1 + 1 

♦ JGZ , -1 (N , 10) 

! N = NMODES 

! I = 1 

TABLE ( NI = " NMODES " , N J=" NMODES " ) : XTKX 
♦LABEL 20 

! K = DS ,"I" ,1 , lP'VLIB" VIBR EVAL 1 1) 

I*"I" : J="I" : "K" 

! 1 = 1 + 1 
♦JGZ , -1 (N , 20) 

♦ IFC’IDSP" NE 0): ♦GOTO 85 
0UTL=22 : INLI=22 

DX = PROD ( DAMP, X) 

"QLIB" XTDX = XTY(X.DX) 

♦LABEL 85 

♦IF("IDMD" NE 0): ♦RETURN 
! N = NMODES : ! I = 1 

♦IF("IDSP" NE 0): TABLE (NI=" NMODES" , NJ= "NMODES") : "QLIB" XTDX 


★ IFC'IDSP" EQ 0): TABLEAU : "QLIB" XTDX 
★LABEL 90 

! D = DS , "I" , 1 , 1 ( M VLIB" DMPD MASK MASK MASK) 

I-" I" : J="I» : "D" 

! I = I + 1 
★JGZ, -1 (N ,90) 

★RETURN 
★LABEL 100 
$ 

$ FULL REDUCED MATRICES (X NOT EIGENVECTORS) 

$ 

★ IFC'DRFO" NE RITZ): ★GOTO 200 
0UTLIB=22 : INLIB=22 

DEFINE K - 1 K SPAR 
DEFINE M = 1 "MNAME" 

KX = PROD(K ,X) 

MX = PRQD(M , X) 

"QLIB" XTKX = XTY(X.KX) 

"QLIB" XTMX = XTY (X , MX) 

★ IFC'IDSP" NE 0): ★GOTO 130 
0UTL=22 : INLI=22 

DX = PROD(DAMP,X) 

"QLIB" XTDX * XTY (X , DX) 

★LABEL 130 

★ IFC'IDMD" NE 0): ★RETURN 
! N = NMODES : ! I = 1 

★ IFC’IDSP" NE 0): TABLE(NI="NMODES" , NJ="NMODES") : ’’QLIB' 1 XTDX 

★ IFC'IDSP" EQ 0): TABLE, U : "QLIB" XTDX 
★LABEL 140 

! D * DS , "I" , 1 , 1 ("VLIB" DMPD MASK MASK MASK) 

I="I" : J= ,, I" : "D" 

! I = I + 1 
★ JGZ, -1 (N , 140) 

★RETURN 
★LABEL 200 
$ 

$ SIMPLE DIAGONAL CASE (X EIGENVECTORS) 

$ 

TABLE (NI=1 ,NJ=" NMODES") : XTKX : TRAN (SOUR=E) 

TABLE(NI=1 ,N J=" NMODES") : XTMX : J=1,"NM0DES" : 1.0 

★ IFC'IDSP" NE 0): ★GOTO 210 
0UTL=22 : INLI*22 

DX = PROD (DAMP ,X) 

"QLIB" XTDX = XTYD (X , DX) 

★LABEL 210 

★ IFC'IDMD" NE 0): ★RETURN 
! N = NMODES : ! I = 1 

★IFC'IDSP" NE 0): TABLE(NI=1 , NJ="NMODES") : "QLIB" XTDX 
★IFC'IDSP" EQ 0): TABLE, U : "QLIB" XTDX 
★LABEL 220 

! D = DS , "I " , 1 , 1 ("VLIB" DMPD MASK MASK MASK) 

1=1 : J="I" : "D" 

! 1 = 1 + 1 
★JGZ, -1 (N , 220) 

★RETURN 

★END 

Runstream TR PLOT 

TR PLOT is a utility runstream for producing plots of response quantities as a function of 
time. Its use is demonstrated in the stepped beam example runstream. 

$ 

$ (TR PLOT) 

$ 

★xqT ui 
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$ PLOTS TRANSIENT, TIME HISTORY DATA PRODUCED BY DR/TRI 

$ 

♦REGISTER EXCEPTIONS TLIB 
♦REGISTER STORE ("TLIB" TR REG 1 1) 

♦REGISTER RETRIEVE! "TLIB" TR REG 1 1) 

$ DEFAULT REGISTER ASSIGNMENTS 


! INLIB=1 
I IDJK = ’ NONE 


! IDQ = ’NONE 
! YNAME = ’DISP 


! N3 = 1 

! TITLE = 'TITLE 
! ID = 1 
! OPT=0 

♦DATA , OPT(TRPLOT OPTIONS) 

$ 

■ NSl = TOC.NIC'INLIB" HIST "YNAME"" N3” MASK) 

! NW1 = TOC.NWDSCINLIB" HIST "YNAME" "N3" MASK) 
i NBLK = TOC , NBLOCKS ( " INLIB" HIST "YNAME" "N3" MASK) 

! NJBL = TOC ,NJ("INLIB" HIST "YNAME" "N3" MASK) 
i NSl ; ! NW1 : ! NBLK : 1 NJBL 
! NSTE=NW1/NS1 
! NPPT=NSTE 

, DT=DS,1 , 1 , 1(" INLIB” DT MASK MASK MASK) $ TIME STEP 
» TIDQ = TOC , IERR(" INLIB" SEL "IDQ” MASK MASK) 
i TIDJ = TOC, IERR("INLIB" SEL ”IDJK” MASK MASK) 
i tTIT = TOC, IERRC "INLIB” "TITLE” MASK MASK MASK) 
♦0NLINE=0 
♦XQT AUS 

TABLE(NI=1 ,NJ="NPPT") : "TLIB” XTAB 
DDATA="DT" 

J=1 , "NPPT” : 0.0 

$ 

$ LOOP OVER ALL RESPONSE QUANTITIES 

I nJLS = 1-NBLK^NJBL + NSTE $ NJ OF LAST BLOCK 
! KBLK = NBLK - 1 
! DBLS = KBLK+NJBL 
! JBLK = KBLK 
! NSM1 = NSl - 1 
! SBASE = 0 


! NJLS 

; H=1 : * N1-NS1 
♦LABEL 

DEFINE Y = "INLIB" HIST "YNAME” ”N3" 1 1,”JBLK 
TABLE (NJ= "NSTE") : "TLIB" YTAB AUS Tl" 

♦ JZ(KBLK ,L2) 

TRAN(S0UR=Y,SBAS="SBAS",SSKIP="NSM1",ILIM=1,JLIM-"NJBL ) 

♦LABEL L2 

DEFINE Y = "INLIB" HIST "YNAME" "N3" 1 "NBLK , NBLK 
tart f it ; "TLIB" YTAB AUS II 

TRAN(S0UR=Y,SBAS="SBAS",SSKIP="NSM1",ILIM=1,JLIM="NJLS",DBAS="DBLS") 

! SBASE = SBASE + 1 
1 11 = 11+1 
♦JGZ, -1 (N1 ,L1) 

$ GENERATE AN X,Y PLOT FOR EACH RESPONSE QUANTITY 


$ 

♦0NLINE=1 
♦XQT PXY 

RESET DEVICE=META 
RESET NDEV=4014 
FONT XNUM=1 : FONT YNUM=1 
FONT XLAB=1 : FONT YLAB=1 
FONT TEXT= 1 
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X - "TLIB" XTAB 
XLABEL 1 TIME (SECONDS) 

♦ JNZ(TIDJ, 110) 

YLF0RMAT,72> 

; (4H J =, 12, 9H JOINT = ,I5,8H COMP = ,12, 8H HIST = ,A4,6H ID - 16) 

♦LABEL 110 ' J 

♦ JNZ(TIDQ , 120) 

YLFORMAT, 72> 

J (4H J” , 12, IX, A4,6H GRP= ,I2,6H IND= ,I5,7H COMP= ,A4,5H ID* 16) 
♦LABEL 120 

XAXIS=3, 5 , 10 
YAXIS=4 ,5,10 
TP0S=O , 0 

! 11*1 : ! N1=NS1 
*JNZ(TTIT,L3) 

TEXT = "TITLE' 1 

♦LABEL L3 

ADVANCE 

BOUNDARIES * .01 .99 .04 .1 
* JNZ(TTIT,L4) 

PLOT TEXT 

♦LABEL L4 

BOUNDARIES* . 0 1 .99 .15 .85 
Y = "TLIB" YTAB AUS "II" 

*JNZ(TIDJ, 210) 

' J0INT = DS,1,"I1'M(»INLIB» SEL "IDJK" MASK MASK) 

! COMP = DS, 2, "II" , 1 ("INLIB" SEL "IDJK" MASK MASK) 

YLABEL "II" "JOINT" "COMP" "YNAME" "ID" 

♦LABEL 210 

♦JNZ(TIDQ,220) 

' ENAHE * DS , 1 , " 1 1 " , 1 ( " INLIB " SEL "IDQ" MASK MASK) 

! EGRP = DS,2, "I 1",1 ("INLIB" SEL "IDQ" MASK MASK) 

! EINDX ~ DS,3, "II" , 1( "INLIB" SEL "IDQ" MASK MASK) 

? EC0MP * DS,4, "II" , 1 ("INLIB" SEL "IDQ" MASK MASK) 

YLABEL "II" "ENAME" "EGRP" "EINDX" "ECOMP" "ID" 

♦LABEL 220 
INIT 

PLOT CURV 
! 11 * 11+1 
♦ JGZ , -1 (N1 , L3) 

♦XQT U1 

♦REGISTER RETRIEVE ("TLIB" TR REG 1 1) 

♦FREE "TLIB" 

♦RETURN 

♦END 


Runstream TR VECTORS 


runslream' l™ ™ VECT0 “ e enerat “ bas « “Mors by calling the system eigensolver, calling 
unstream TR RITZ, using the static mode method, or by other experimental techniques. 

$ — 

$ (TR, VECTORS) 

$ 

♦XQT AUS 
R = RIGID(l) 


COMPUTE VECTORS FOR USE IN DYNAMIC ANALYSIS 


CR = PROD ( "MNAME" , R) 

Z = NDDF, 1 (CR) 

• NDDF = DS, 1,1, 1(1 z AUS 1 1) 
! NDDF 


♦IF ("METHOD" NE MODE) : *GOTO 100 
♦XQT E4 

RESET NMODES* " NMODES " 

RESET M*" MNAME" 

RESET NDDF* "NDDF" 

RESET CONV*"CONV" 
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♦DCALL,0PT(E4 PARAMETERS) 

♦GOTO 200 

$ 

♦LABEL 100 

♦IF ("METHOD” NE RITZ) : ♦GOTO 105 
! MAXHZ = 1.03E+10 
! NMD = NMODES 
! SCALE =1.0 
! AFLIB = 1 
! VLIB = 1 
♦DCALL(TR RITZ ) 

♦GOTO 200 
$ 

♦LABEL 105 

♦IF ( "METHOD" NE OLD): ♦GOTO 110 
♦XQT DCU 

COPY 3 1 VIBR MODE 
COPY 3 1 VIBR EVAL 
♦GOTO 200 
♦LABEL 110 

♦IFC'METHOD” NE ONES): ♦GOTO 115 
♦DCALLCTEST NEB3) 

♦GOTO 200 
♦LABEL 115 

♦IFC'METHOD” NE STAT) : ♦GOTO 120 
♦XQT E4 

RESET NMODES= "NMODES" 

RESET M="MNAME" 

RESET NDDF=”NDDF" 

RESET CONV="CONV" 

♦DCALL,0PT(E4 PARAMETERS) 

♦XQT DRSI 
♦XQT SSOL 

$ ORTHOGONALIZE STATIC SOLUTION AND APPEND TO SET OF MODE SHAPES 
♦DCALL(TR.GRAM) 

$ MAKE THE VECTORS ORTHOGONAL WITH RESPECT TO BOTH K AND M 
*DCALL(TR DIAG) 

♦XQT VPRT 
PRINT SI AUS 
♦XQT DCU 
TOC 1 

PRINT 1 VIBR EVAL 
♦LABEL 120 

♦IFC'METHOD” NE UMOT) : *GQ TO 200 
♦XQT AUS 

Z = NDDF ,1,2 (CR) 

♦XQT DRSI 
RESET CDN=2 
♦XQT AUS 
UDF = 1 
SSPREP (K , 2) 

♦XQT DCU 

CHANGE 1 BNF MASK MASK MASK VIBR MODE 1 1 
$ UNCOUPLE THE SYSTEM 
♦DCALL(TR.DIAG) 

$ END OF METHODS OPTIONS 

♦LABEL 200 

♦END 


Runstream TR GRAM 

TR GRAM is a utility runstream for performing a Gram-Schmidt orthogonalization of a set 
of vectors. 

$ 

$ (TR,GRAM) - PERFORM GRAM-SCHMIDT PROCESS TO M-QRTHOGDNALIZE 
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$ STAT DISP WITH RESPECT TO VIBR MODE AND THEN 

$ REPLACE THE LAST VIBR MODE WITH THE NEW VECTOR 

$ 

★XQT AUS 

! NMM1 55 NMODES - 1 
DEFINE X = VIBR MODE 1 1 1,"NMM1" 

DEFINE S * STAT DISP 1 1 
DEFINE M = "MNAME" 

INLIB=10 : QUTLIB=10 

$ ORTHOGONALIZE THE STATIC SOLUTION WITH RESPECT TO THE MODE SHAPES 
MX = PRQD(M, X) 

XTMS = XTY (MX ,S) 

A = CBR(X, XTMS) 

SI = SUM(S, -1.0 A) 

$ NOW SCALE THE VECTOR 
MS = PROD (M , SI) 

STMS * XTY(Sl.MS) 

! STMS = DS, 1,1, 1(10 STMS AUS 11):! STMS = STMS**. 5 : ! STMS * 1.0/STMS 
TEMP MODE 1 1 = UN I ON ("STMS " S1,X) 

*XQT DCU 

CHANGE 10 TEMP MODE 1 1 VIBR MODE 1 1 
COPY 10 1 VIBR MODE 1 1 
★DELETE 10 
♦END 

Runstream SENS DVUP 

SENS DVUP is a utility runstream for updating the design variable registers based on the 
data sets X ADS and XNAME ADS. It is always called immediately before calling MODEL so the 
current values of the design variables are available for use. 

$ 

$ (SENS DVUP) - UPDATE DESIGN VARIABLE REGISTERS FROM DATASET 

$ 

*XQT U1 

! N = TOC,NJ (1 XNAME ADS 1 1) 

’ 1*1 
★LABEL 10 

! RNAME = DS,1,"I",1(1 XNAME ADS 1 1) 

! RVAL * DS, 1 , "I" , 1(1 X ADS 1 1) 

! "RNAME" = "RVAL" 

! "RNAME" 

! 1 = 1 + 1 
* JGZ, -1 (N , 10) 

★ END 

Runstream TR DXDV 1 

The TR DXDV n runstreams implement the different sensitivity methods. The structure of 
all these runstreams is similar. In each case there is a loop over the designated design variables, 
and sensitivities are calculated of the required response quantities at the set of critical points. 
Within this loop there is at least one call to runstream MODEL to form a perturbed design, 
a call to form a set of new reduced equations, and a call to processor DRX to integrate the 
reduced equations in time. Runstream TR DXDV 1 implements the forward difference method 
with either fixed or updated basis vectors. 

$ 

$ (TR DXDV 1) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
$ USING THE FORWARD DIFFERENCE OPERATOR AND 

$ EITHER FIXED OR UPDATED MODES 

$ UPDATE HISTORY 

$ 6/28/88 WHG - MODIFIED FOR VELD, ACCE, STRESSES 

$ 

★XQT U1 

♦REGISTER ST0RE(1 DXDV REGISTERS 1 1) 
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♦REGISTER RETRIEVE ( 1 DXDV REGISTERS 1 1) 

♦RGI 

FDCH = .001 
FDMCH = .0001 
XLIB = 5 
OPT = 0 
RUB 3 14 
DRMETH0D=0 
DXMD-UPDATED 

♦DCALL , OPT (DXDV PARAMETERS) 

♦SHOW 

♦DCALL (TR , DPREP ) 

$ 

$ LOOP OVER ALL DESIGN VARIABLES 
$ 

! NDV * TOC ,NJ ( 1 X ADS 1 1) 

! NCNT = NDV 
! IDV - 1 
♦LABEL 10 
♦xqT U1 

! IFLG = DS , 1 , "IDV” ,1(1 XFLG ADS 1 1) 

♦JZ(IFLG , 100) 

♦LIBS "XLIB” 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦XQT DCU 

COPY "XLIB” 1 XNAME ADS 1 1 
♦XQT U3 
RP2 

FORMAT 1 ’ (1H1 , 20X , 41HBEGINNING SENSITIVITY CALCULATION. DV NO., 13) 
PRINT Cl) "IDV" 

♦XQT AUS 

DEFINE X * "XLIB" X ADS 1 1 
TABLE(NI=1,NJ 3 "NDV") : X ADS 1 1 
J=1 , "NDV" : 1.0 
J="IDV" : "FDCH” 

TRAN (SOURCES, OPERATION-MULT) 

$ CHECK FOR TOO SMALL A STEP 
! X = DS , 1 , "IDV” , 1 ( "XLIB" X ADS 1 1) 

! DX = FDCH^X 

♦IF(”DX” GT "FDCHM”) : *GOTO 20 
! DX = FDMCH 
! X = X + FDMCH 
TABLE, U : X ADS 1 1 
OPER = XSUM 
J="IDV" : "X" 

♦LABEL 20 
♦ CALL ( SENS ,DVUP) 

$ FORM PERTURBED MODEL 

♦0NLINE=0 

♦DCALL (MODEL) 

♦ONLINE 3 1 
♦XQT DCU 


COPY 

"XLIB" 

1 

TIME 

COPY 

"XLIB" 

1 

CA 

COPY 

"XLIB" 

1 

DMPD 


♦IF(”DXMD" NE FIXE): *G0T0 30 
COPY "XLIB" 1 VIBR MODE 
♦DCALL (TR.DIAG) 

♦GO TO 40 
♦LABEL 30 

♦IF("DXMD" NE UPDA): *0010 40 
♦DCALL (TR, VECTORS) 

♦LABEL 40 
♦XQT AUS 

DEFINE X = VIBR MODE 1 1 1,"NM0DES" 
DEFINE F = APPL FORC 1 
XTF = XTY(X.F) 



♦DCALL (TR , REDM) VLIB=1 
♦XQT DRX 

DTEX (DT= " DT " , METHOD= " DRMETHOD 11 , NTERMS= " NTERMS " ) 

TR1 (QXLIB=1 , QX1L=1 , QX2L=1 , T2=' , T2" , LB="BLKSIZE") 

♦DCALL (TR,DBACK 1) 

$ 

$ COMPUTE DERIVATIVES USING FORWARD DIFFERENCE OPERATOR 
$ 

! OVDX = l.O/DX 
! MOVD = - OVDX 
•1=1: ! N = NBCK 

♦LIBS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦LABEL 80 
♦XQT AUS 

• NM = DS , "I" , 1 , 1(1 BACK LIST 1 1) 

! I ERR = TOC , IERR(1 SEL "NM" MASK MASK) 

♦ JNZ (IERR, 90) 

DEFINE CPI - "XLIB" CRPT "NM" 

DEFINE CPO = CRPT "NM" 

DXDV "NM" "IDV" * SUM ("OVDX " CPI "MOVD" CPO) 

♦XQT DCU 

PRINT 1 DXDV "NM" "IDV" 

♦LABEL 90 
! I = I + 1 
♦ JGZ , -1 (N , 80) 

ERASE "XLIB" 

♦LABEL 100 
! IDV = IDV + 1 
♦JGZ , -1 (NCNT, 10) 

$ 

♦XQT U1 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦RETURN 

♦END 

Runstream TR DXDV 3 


Runstream TR DXDV 3 implements the fixed-mode semianalytical sensitivity method De- 
pending on the call to runstream TR DBACK, either the mode displacement or mode acceleration 
method is used to recover the physical sensitivities. 

$ 

$ (TR DXDV 3) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
$ SEMI ANALYTICALLY 

$ 

♦XQT U1 #■ 

♦REGISTER STORE (1 DXDV REGISTERS 1 1) 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦RGI 

FDCH = .001 
FDMCH = .0001 
XLIB = 5 
OPT = 0 
RLIB = 14 
DRMETH0D=0 

♦DCALL, OPT (DXDV PARAMETERS) 

♦SHOW 

$ 

$ INITIALIZATION FOR DERIVATIVE CALCULATIONS 
$ 

♦DCALL (TR,DPREP) 

♦XQT U1 
$ 

$ LOOP OVER ALL DESIGN VARIABLES 
$ 

! NDV = T0C,NJ(1 X ADS 1 1) 
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! NCNT * NDV 
! IDV = 1 
♦LABEL 10 

♦LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦XQT U1 

! IFLG = DS , 1 , "IDV” , 1 ( "XLIB" XFLG ADS 1 1) 

♦JZ(IFLG, 100) 

♦XQT DCU 

COPY "XLIB” 1 XNAME ADS 1 1 
♦XQT U3 
RP2 

FORMAT 1 * ( 1H1 , 20X , 4 1HBEG INNING SENSITIVITY CALCULATION. DV NO. ,13) 
PRINT(l) "IDV” 

♦XQT AUS 

DEFINE X = "XLIB" X ADS 1 1 
TABLE(NI=1,NJ="NDV") : X ADS 1 1 
J=1 , "NDV" : 1.0 
J="IDV" : "FDCH" 

TRAN (SOURCE=X , OPERAT I ON=MULT ) 

$ CHECK FOR TOO SMALL A STEP 
! X = DS , 1 , "IDV" , 1 ("XLIB" X ADS 1 1) 

! DX = FDCH^X 

♦IF("DX" GT "FDCHM") : *GOTO 20 
! DX = FDMCH 
! X = X + FDMCH 
TABLE ,U : X ADS 1 1 
OPER = XSUM 
J-"IDV" : "X" 

♦LABEL 20 
♦CALL (SENS, DVUP) 

$ FORM PERTURBED MODEL 

♦ONLINE^O 

♦DCALL (MODEL) 

♦ONLINE- 1 
♦XQT AUS 

DEFINE X = "XLIB" VIBR MODE 1 1 1,"NM0DES" 

DEFINE KO = "XLIB" K SPAR 

DEFINE Ki = K SPAR 

DEFINE MO = "XLIB" "MNAME" 

DEFINE Ml = "MNAME" 

DEFINE DO = "XLIB" DAMP SPAR 

DEFINE D1 = DAMP SPAR 

DEFINE FO = "XLIB" APPL FORC 1 

DEFINE FI = APPL FORC 1 

! IDSP = TOC , IERR(1 DAMP SPAR MASK MASK) 

! OVDX - 1.0/DX 
! MOVD = -OVDX 

DKDV = SUM ("OVDX" Kl "MOVD" KO) 

DMDV = SUM ("OVDX" Ml "MOVD" MO) 

DFDV = SUM("DVDX" FI "MOVD" FO) 

♦ IFCIDSP" Eq 0): DDDV = SUM("OVDX" D1 "MOVD" DO) 

DKX = PROD (DKDV , X) 

DMX = PROD (DMDV, X) 

♦ IFCIDSP" EQ 0): DDX = PROD(DDDV,X) 

RDKX - XTY (X ,DKX) 

RDMX = XTY (X ,DMX) 

♦ IFCIDSP" EQ 0): RDDX = XTY(X,DDX) 

XTF AUS = XTY (X, DFDV) 

♦XQT DCU 

COPY "XLIB" 1 TIME 

COPY "XLIB" 1 CA 

COPY "XLIB" 1 DT AUS 

COPY "XLIB" 1 DTEX AUS 

COPY "XLIB" 1 DCON AUS $ CONSTANTS FOR NEWMARK METHOD 

$ 

$ FORM THE RIGHT-HAND-SIDE PSEUDO LOAD VECTOR 
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$ 

♦XQT DRX 

BACK (LRZ="BLKSI ZE" ,PRINT=0) 

T - -1.0 RDMX : Y = "XLIB" QX2 AUS 
♦IF ( "IDSP" EQ 0): T » -1.0 RDDX : Y = "XLIB" QX1 AUS 
T == -1.0 RDKX : Y = "XLIB" QX AUS 
Z = FH AUS 
♦XQT DRX 

TR1 (QXLIB-1 , QX1LIB=1 , QX2LI3=1 ,T2= M T2" , FHLIB=1 , LB="BLKSIZE" ) 

$ 

$ BACK TRANSFORM FOR NECESSARY SENSITIVITIES OF PHYSICAL 
$ QUANTITIES 
$ 

♦DCALL (TR , DBACK ,4) 

♦XQTC DCU 
ERASE 1 
♦LABEL 100 
! IDV * IDV + 1 
♦JGZ ,-l (NCNT, 10) 

$ 

♦LIBS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦XQT U1 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦RETURN 

♦END 

Runstream TR DXDV 5 

Runstream TR DXDV 5 implements the overall central difference method using either fixed 
or updated basis vectors. 

$ 

$ (TR DXDV 5) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
$ USING TWO POINT CENTRAL DIFFERENCE OPERATOR 

$ WITH UPDATED OR FIXED MODES 

$ 

♦XQT U1 

♦REGISTER STDRE ( 1 DXDV REGISTERS 1 1) 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦ RGI 

FDCH = .001 
FDMCH = .0001 
XLIB = 5 
YLIB = 6 
OPT = 0 
RLIB = 14 
DRMETH0D=0 
DXMD=UPDATED 

♦DCALL, OPT (DXDV PARAMETERS) 

♦SHOW 

$ 

$ INITIALIZATION FOR DERIVATIVE CALCULATIONS 

$ 

♦DCALL (TR , DPREP ) 

$ 

$ LOOP OVER ALL DESIGN VARIABLES 

$ 

! NDV = T0C,NJ(1 X ADS 1 1) 

! NCNT = NDV 
! IDV = 1 
♦LABEL 10 
♦XQT U1 

! IFLG = DS , 1 > " IDV" ,1(1 XFLG ADS 1 1) 

*JZ(IFLG , 100) 

♦XQT U3 
RP2 
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FORMAT V (1H1.20X.41HBEGINNING SENSITIVITY CALCULATION. DV NO., 13) 
PRINT(l) "IDV" 

♦LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦XQT AUS 
! IS * 1 
! NST = 2 
! SIGN =1.0 
$ 

$ DO ANALYSIS FOR BOTH POSITIVE AND NEGATIVE STEPS 
$ 

♦LABEL 15 
♦XQT DCU 

COPY "XLIB" 1 XNAME ADS 1 1 
♦XQT AUS 

DEFINE X = "XLIB" X ADS 1 1 
TABLE(NI=1 ,NJ= ,, NDV") : X ADS 1 1 
TRAN (SOURCE=X) 

$ DIFFERENCE APPROPRIATE DESIGN VARIABLE 
! X = DS , 1 , "IDV" , 1 ("XLIB" X ADS 1 1) 

! DX = FDCH*X 

♦IF("DX" GT "FDCHM") : ♦GOTO 20 
! DX = FDMCH 
♦LABEL 20 
! X = DX+SIGN + X 
TABLE ,U : X ADS 1 1 
OPER = XSUM 
J="IDV" : "X" 

♦CALL ( SENS , DVUP ) 

$ FORM PERTURBED MODEL 

♦0NLINE=O 

♦DCALL (MODEL) 

♦0NLINE=1 
♦XQT DCU 

COPY "XLIB" 1 TIME 
COPY "XLIB" 1 CA 
COPY "XLIB" 1 DMPD 
♦IF( "DXMD" NE FIXE): ♦GOTO 30 
COPY "XLIB" 1 VIBR MODE 
♦DCALL (TR,DI AG) 

♦GO TO 40 
♦LABEL 30 

♦IF("DXMD" NE UPDA): *GOTO 40 
♦DCALL (TR, VECTORS) 

♦LABEL 40 
♦XQT AUS 

DEFINE X = VIBR MODE 1 1 1,"NM0DES" 

DEFINE F = APPL FORC 1 
XTF = XTY(X.F) 

♦DCALL (TR.REDM) VLIB=1 
♦XQT DRX 

DTEX(DT="DT" , METHOD="DRMETHDD" ,NTERMS="NTERMS" ) 

TR1 (QXLIB=1 , T2="T2" , QX1LIB=1 , QX2LIB=1 , LB="BLKSIZE" ) 

♦DCALL (TR.DBACK 1) 

! SIGN = -1.0 

♦LIBS "YLIB" 2 3 4 1 5 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
* JGZ , -1 (NST , 15) 

$ 

$ COMPUTE DERIVATIVES USING CENTRAL DIFFERENCE OPERATOR 

$ 

‘ TWDX = 2 . O^DX 
! OVDX = 1.0/TWDX 
! MOVD = - OVDX 
! I = 1 : ! N = NBCK 

♦LIBS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦LABEL 80 
♦XQT AUS 
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■ NM = DS, "I" , 1,1(1 back LIST 1 1) 

! IERR = TOC , IERR( 1 SEL "NM" MASK MASK) 

♦JNZ(IERR, 90) 

DEFINE CPI = "XLIB" CRPT "NM" 

DEFINE CPO = "YLIB" CRPT "NM" 

DXDV "NM" "IDV" = SUM ( "OVDX" CPI "MOVD" CPO) 

♦xqT DCU 

PRINT 1 DXDV "NM" "IDV" 

♦LABEL 90 
! I = I + 1 
♦JGZ,-1(N,80) 

ERASE "XLIB" 

ERASE "YLIB" 

♦LABEL 100 
! IDV = IDV + 1 
* JGZ , -1 (NCNT , 10) 

$ 

♦XQT U1 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦RETURN 

♦END 

Runstream TR DXDV 6 

Runstream TR DXDV 6 implements the semianalytical method with nonzero d^/dx. The 
called procedure TR DPHI determines how the basis vector derivatives are calculated. 

$ 

$ (TR DXDV 6) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
$ SEMI ANALYTICALLY BUT WITH THE EFFECT OF CHANGING 

$ MODES INCLUDED 

$ 

♦XQT U1 

♦REGISTER STORE (1 DXDV REGISTERS 1 1) 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦RGI 

FDCH = .001 
FDMCH = .0001 
XLIB * 5 
DPT = 0 
RLIB = 14 
DRMETH0D=0 

♦DCALL, DPT (DXDV PARAMETERS) 

♦SHOW 

$ 

$ INITIALIZATION FOR DERIVATIVE CALCULATIONS 

$ 

♦DCALL ( TR , DPREP) 

♦XQT U1 

$ 

$ LOOP OVER ALL DESIGN VARIABLES 

$ 

! NDV * TOC ,NJ ( 1 X ADS 1 1) 

! NCNT = NDV 
! IDV - 1 
♦LABEL 10 

♦LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 IT 18 19 20 
♦XQT U1 

! IFLG = DS , 1 , "IDV" , 1 ("XLIB" XFLG ADS 1 1) 

♦JZ(IFLG.IOO) 

♦XQT DCU 

COPY "XLIB" 1 XNAME ADS 1 1 
♦XQT U3 
RP2 

FORMAT 1 ’ ( 1H1 , 20X , 41HBEGINNING SENSITIVITY CALCULATION. DV NO 13) 

PRINT ( 1) "IDV" 
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♦XQT AUS 

DEFINE X = "XLIB" X ADS 1 1 
TABLE(NI=1 ,NJ="NDV") : X ADS 1 1 
J=1 , "NDV" : 1.0 
J="IDV" : "FDCH” 

TRAN (S0URCE=X , 0PERATI0N=KULT) 

$ CHECK FOR TOO SMALL A STEP 
! X = DS , 1 , "IDV" , 1 ("XLIB" X ADS 1 1) 

! DX = FDCH*X 

♦ IFO'DX" GT "FDCHM" ) : ♦GOTO 20 
! DX = FDMCH 
! X = X + FDMCH 
TABLE, U : X ADS 1 1 
OPER = XSUM 
J="IDV" : "X" 

♦LABEL 20 
♦CALL (SENS, DVUP) 

! OVDX = 1.0/DX 
! MOVD = -OVDX 
$ FORM PERTURBED MODEL 
♦QNLINE=0 
♦DCALL (MODEL) 

$ CALCULATE DERIVATIVES OF MODES SHAPES 
♦XQT AUS 

DEFINE KO = "XLIB" K SPAR 

DEFINE K1 = K SPAR 

DEFINE MO = "XLIB" "MNAME" 

DEFINE Ml = "MNAME" 

DKDV = SUM("OVDX" K1 "MOVD" KO) 

DMDV = SUM ( "OVDX" Ml "MOVD" MO) 

♦DCALL (TR,DPHI , 3) 

♦XqT DCU 

COPY "XLIB" 1 DT AUS 
COPY "XLIB" 1 DTEX AUS 

COPY "XLIB" 1 DCON AUS $ CONSTANTS FOR NEWMARK METHOD 
♦XQT AUS 

DEFINE XO = "XLIB" VIBR MODE 1 1 1,"NM0DES" 

DEFINE KO = "XLIB" K SPAR 

DEFINE K1 = K SPAR 

DEFINE MO = "XLIB" "MNAME" 

DEFINE Ml = "MNAME" 

DEFINE DO = "XLIB" DAMP SPAR 

DEFINE D1 = DAMP SPAR 

DEFINE FO = "XLIB" APPL FORC 1 1 

DEFINE FI = APPL FORC 1 1 

DEFINE DXDV = DXDV AUS "IDV" 

$ 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE STIFFNESS MATRIX 

$ 

DKX1 = PROD ( DKDV, XO) 

DKX2 = PROD CKO ,X0) 

DKX3 = PROD (KO, DXDV) 

XDK1 = XTY(XO ,DKX1) 

XDK2 = XTY (DXDV ,DKX2) 

XDK3 = XTY (XO ,DKX3) 

XTMP = SUM(XDK1 , XDK2) 

XDKX = SUM (XTMP , XDK3) 

$ 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE MASS MATRIX 

$ 

DMX1 = PROD ( DMDV, XO) 

DMX2 = PROD (MO , XO) 

DMX3 = PROD (MO, DXDV) 

XDM1 = XTY (XO ,DMX1 ) 

XDM2 = XTY (DXDV ,DMX2) 

XDM3 = XTY (XO , DMX3) 


XTMP = SUM (XDM 1 , XDM2 ) 

XDMX = SUM ( XTMP, XDM3) 

$ 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE DAMPING MATRIX 

$ 

• IDSP * TOC , IERR( 1 DAMP SPAR MASK MASK) 

♦IF ( "IDSP" NE 0): ♦GOTO 30 
DDDV * SUMO'OVDX" D1 "MOVD" DO) 

DDX1 = PROD (DDDV, XO) 

DDX2 - PROD (DO, XO) 

DDX3 = PROD(DO,DXDV) 

XDD1 = XTY (XO.DDXl) 

XDD2 =* XTY (DXDV ,DDX2) 

XDD3 * XTY(X0,DDX3) 

XTMP =* SUM(XDD1 ,XDD2) 

XDDX = SUM (XTMP ,XDD3) 

♦LABEL 30 

$ 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE FORCE VECTOR 

$ 

DFDV * SUMC'OVDX" FI "MOVD" FO) 

XF1 = XTY (DXDV, FO) 

XF2 * XTY(XO, DFDV) 

XTF AUS * SUM(XF1 , XF2) 

♦XQT DCU 
PRINT 1 XTF AUS 
COPY "XLIB" 1 TIME AUS 
CDPY "XLIB" 1 CA AUS 

$ 

$ FORM THE RIGHT-HAND-SIDE PSEUDO LOAD VECTOR 

$ 

♦XQT DRX 

BACK (LRZ= s "BLKSIZE" ,PRINT=0) 

T - -1.0 XDMX : Y = "XLIB" QX2 AUS 
♦IF("IDSP" EQ 0): T * -1.0 XDDX : Y * "XLIB" QX1 AUS 
T = -1.0 XDKX : Y = "XLIB" QX AUS 
Z = FH AUS 
♦XQT DRX 

TR1 (QXLIB*1 , QX1LIB-1 , QX2LIB=1 , T2="T2" , FHLIB^l , LB“"BLKSIZE" ) 
♦XQT DCU 
TOC 1 

$ 

$ BACK TRANSFORM FOR NECESSARY SENSITIVITIES OF PHYSICAL 
$ QUANTITIES 
$ 

♦DCALL(TR,DBACK , 3) 

♦XQTC DCU 
ERASE 1 
♦LABEL 100 
! IDV = IDV + 1 
♦JGZ,-1(NCNT, 10) 

$ 

♦LIBS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 

♦0NLINE*1 

♦XQT U1 

♦REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

♦RETURN 

♦END 

Runstream TR DPREP 


TR DPREP is a utility runstream used by all the sensitivity calculation runstreams. 
task is to locate the critical points for all required response quantities. 

$ 

$ (TR, DPREP) - PREPARATION FOR SENSITIVITY CALCULATIONS 


Its main 
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$ 

$ 

$ FORM CRITICAL POINT TABLES FOR RESPONSE QUANTITIES 
$ 

! NBCK = TOC ,NI ("QLIB" BACK LIST 1 1) 

! I = 1 : ! N = NBCK 
■►LABEL 10 

! NM = DS, "I" , 1 , 1 ("QLIB" BACK LIST 1 1) 

! IERR = TOC, IERR( "QLIB" SEL "NM" MASK MASK) 
*JNZ(IERR»20) 

♦XQT U10 

CRIT(Y="QLIB" HIST "NM" ,DT="DT" , NCRIT* " NCRIT " , k 

CRPT="QLIB" CRPT "NM" , CRTI="QLIB" CRTI "NM" , k 
PCH= . 25) 

♦XQT DCU 

PRINT "QLIB" CRTI "NM" 

PRINT "QLIB" CRPT "NM" 

♦LABEL 20 
! 1 = 1 + 1 
♦ JGZ , — 1 ( N , 10) 

♦XQT U1 

♦ (E4 PARAMETERS) EOFX 
$ RESET NFCT="NMODES" , NLIM="NMODES" 

RESET NIF="NMODES" 

IFSOURCE= "XLIB" VIBR MODE 1 1 
♦EOFX 
♦END 


Runstream TR DBACK 1 

TR DBACK n runstreams implement the different procedures for recovering the physical 
sensitivities. They all rely heavily on runstream TR CRPT which recovers a specific physical 
quantity at the critical points. Runstream TR DBACK 1 recovers the sensitivities with the mode 
displacement method and is used in the overall finite difference procedures. 



$ (TR DBACK 1) - BACK TRANSFORMATION FOR DERIVATIVES USING MODE 
$ DISPLACEMENT METHOD 

$ 

! I IB = 1 : ! NNB = NBCK 
♦LABEL 10 
♦XQT AUS 

DEFINE X = VIBR MODE 1 1 1,"NM0DES" 

! NM = DS, "IIB" , 1 , 1("XLIB" BACK LIST 1 1) 

! IERR = TOC, IERR( "XLIB" SEL "NM" MASK MASK) 

♦JNZCIERR.200) 

$ 

$ DISPLACEMENTS 
$ 

♦ IFO'NM" NE DISP): ♦GOTO 30 
DEFINE IDJK = "XLIB" SEL DISP 
TMAT VMOD = SVTRAN ( IDJK , X) 

♦DCALL (TR , CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VMOD CNAME=DISP Q=QX 
♦GOTO 200 
♦LABEL 30 
$ 

$ VELOCITIES 

$ 

♦ IFO'NM" NE VELO) : ♦GOTO 50 
DEFINE IDJK = "XLIB" SEL VELO 
TMAT VVEL = SVTRAN ( IDJK , X) 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VVEL CNAME=VELO Q=QX1 
♦GOTO 200 
♦LABEL 50 
$ 
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$ ACCELERATIONS 
$ 

♦ IFO'NM" NE ACCE): ♦GOTO 70 
DEFINE IDJK = "XLIB" SEL ACCE 
TMAT VACC = SVTRAN ( IDJK , X) 

♦DC ALL (TR , CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VACC CNAME=ACCE Q=QX2 
♦GOTO 200 
♦LABEL 70 
$ 

$ REACTIONS 

$ 

♦IFO'NM" NE REAC) : ♦GOTO 90 
♦GOTO 200 
♦LABEL 90 
$ 

$ STRESSES 

$ 

♦ IFO'NM" NE STRE): ♦GOTO 110 
♦XQT ES 

RESET OPER=T 

IDQ = "XLIB" SEL STRESS 

U * VIBR MODE 1 1 1 , "NMODES" 

T = TMAT VSTRE 

♦DCALL ( TR , CRPT ) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VSTRE CNAME-STRE Q=QX 

♦GOTO 200 

♦LABEL 110 

♦LABEL 200 

! IIB = IIB + 1 

♦JGZ, -1 (NNB, 10) 

! IIB = FREEQ : ! NNB = FREE() : ! LIB1 = FREEO : ! LIB2 = FREEO 
! LIB3 = FREEO : ! TNAME = FREEO : ? CNAME = FREEO 
! Q = FREEO 
♦END 

Runstream TR DBACK 2 

Runstream TR DBACK 2 recovers sensitivities in the fixed-mode, mode displacement version 
of the semianalytical method. 

$ 

$ (TR DBACK 2) - BACK TRANSFORMATION FOR DERIVATIVES USING 
$ SEMIANALYTICAL METHOD 

$ UPDATE HISTORY 

$ 

$ 6/22/88 WHG - MODIFIED FOR SEMIANALYTICAL METHOD 
$ 6/21/88 WHG - CREATED FROM ( TR, DBACK, 1) FOR UPDATED MODES 

$ 

! IIB - 1 : ! NNB = NBCK 
♦LABEL 10 
♦XQT AUS 

! NM = DS, "IIB", 1,1 ("XLIB" BACK LIST 1 1) 

! I ERR = TOC, IERRC'XLIB" SEL "NM" MASK MASK) 

*JNZ(IERR,300) 

$ 

$ DISPLACEMENTS 
$ 

♦ IFO'NM" NE DISP): ♦GOTO 30 

♦DCALL (TR , CRPT ) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VMDD CNAME=DISP Q*QX 
♦GOTO 200 
♦LABEL 30 
$ 

$ VELOCITIES 

$ 

♦ IFO'NM" NE VELD): ♦GOTO 50 

♦DCALL (TR, CRPT) LIB1*"XLIB" LIB2="XLIB" LIB3=1 TNAME=VVEL CNAME=VELO Q-QX1 
♦GOTO 200 
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■►LABEL 50 

$ 

$ ACCELERATIONS 
$ 

★IF( "NM" NE ACCE): ★GOTO 70 

★DC ALL ( TR , CRPT ) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VACC CNAME=ACCE Q-QX2 
★GOTO 200 
★LABEL 70 
$ 

$ REACTIONS 
$ 

★IF ("NM" NE REAC) : ★GOTO 90 
★GOTO 200 
★LABEL 90 
$ 

$ STRESSES 

$ 

★IF("NM" NE STRE): ★GOTO 110 
$ FORM [S] [DQ/DV] 

★DCALL(TR,CRPT) LIB1="XLIB" LIB2-"XLIB" LIB3=1 TNAME=VSTRE CNAME=STRE Q-QX 
★XQT DCU 

CHANGE 1 CRPT STRE 1 1 CRPT STR1 1 1 
$ FORM [DS/DV] [Q] 

★XQT ES 
RESET OPER=T 
IDQ = "XLIB" SEL STRESS 
U = "XLIB" VIBR MODES 1 1 1,"NM0DES M 
T = TMAT VSTRE 
★XQT AUS 

DEFINE SO = "XLIB" TMAT VSTRESS 

DEFINE SI = TMAT VSTRESS 

TMAT DSDV = SUMC'OVDX" SI "MOVD" SO) 

★DC ALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSDV CNAME-STRE Q*QX 
★XQT AUS 

DEFINE STR1 = CRPT STR1 

DEFINE STR2 = CRPT STRE 

CRPT STRE - SUM(STR1 , STR2) 

★GOTO 200 
★LABEL 110 
★GOTO 300 
★LABEL 200 
★XQT DCU 

CHANGE 1 CRPT "NM" 1 1 DXDV "NM" "IDV" 1 
COPY 1 "XLIB" DXDV "NM" "IDV" 1 
PRINT "XLIB" DXDV "NM" "IDV" 1 
★LABEL 300 
! IIB = IIB + 1 
★JGZ , -1 (NNB, 10) 

! IIB = FREE ( ) : ! NNB = FREE ( ) : ! LIB1 = FREE() : ! LIB2 = FREE() 

! LIB3 = FREEO : ! TNAME = FREE() : ! CNAME = FREEO 

! q = FREEO 

★END 


Runstream TR DBACK 3 

Runstream TR DBACK 3 recovers sensitivities in the semianalytical method with nonzero 
d<&/dx. 

$ 

$ (TR DBACK 3) - BACK TRANSFORMATION FOR DERIVATIVES USING 
$ SEMIANALYTICAL METHOD WITH CHANGING MODES 

$ UPDATE HISTORY 
<£ 

$ 6/22/88 WHG - MODIFIED FOR SEMIANALYTICAL METHOD 
$ 6/21/88 WHG - CREATED FROM (TR, DBACK , 1) FOR UPDATED MODES 



! I IB = 1 : ! NNB = NBCK 
♦LABEL 10 
♦XQT AUS 

DEFINE DX = DXDV AUS "IDV" 1 1 , "NMODES" 

• NM = DS , "IIB" , 1 , 1 ("XLIB" BACK LIST 1 1) 

! I ERR = TOC , IERRO'XLIB" SEL "NM" MASK MASK) 

♦JNZCIERR.300) 

$ 

$ DISPLACEMENTS 
$ 

♦IF("NM" NE DISP): ♦GOTO 30 
DEFINE IDJK = "XLIB" SEL DISP 
TMAT DVMX = SVTRAN(IDJK ,DX) 

*DCALL(TR,CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VMOD CNAME=DISP Q=QX 
♦XQT DCU 

CHANGE 1 CRPT DISP 1 1 CRPT DSP1 1 1 

♦DCALL (TR , CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DVMX CNAME*DISP Q=QX 
♦XqT AUS 

DEFINE D1 = CRPT DSP1 
DEFINE D2 = CRPT DISP 
CRPT DISP = SUM (D1 , D2) 

♦GOTO 200 
♦LABEL 30 
$ 

$ VELOCITIES 
$ 

♦IF("NM" NE VELO): ♦GOTO 50 
DEFINE IDJK = "XLIB" SEL VELO 
TMAT DVMX = SVTRAN (IDJK , DX) 

♦DCALL (TR » CRPT) LIB1~"XLIB" LIB2="XLIB" LIB3-1 TNAME=VVEL CNAME=VELO Q=QX1 
♦XQT DCU 

CHANGE 1 CRPT VELO 1 1 CRPT VEL1 1 1 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DVMX CNAME=VELO Q=QX1 
♦XQT AUS 

DEFINE VI = CRPT VEL1 
DEFINE V2 = CRPT VELO 
CRPT VELO = SUM(V1,V2) 

♦GOTO 200 
♦LABEL 50 
$ 

$ ACCELERATIONS 
$ 

♦IFO'NM" NE ACCE) : ♦GOTO 70 
DEFINE IDJK = "XLIB" SEL ACCE 
TMAT DVMX = SVTRAN ( IDJK , DX) 

♦DCALL (TR, CRPT) LIB1*"XLIB" LIB2="XLIB" LIB3=*1 TNAME=VACC CNAME=ACCE q=QX2 
♦XQT DCU 

CHANGE 1 CRPT ACCE 1 1 CRPT ACC1 1 1 

♦DCALL (TR , CRPT ) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DVMX CNAME-ACCE Q=QX2 
♦XQT AUS 

DEFINE A1 = CRPT ACC1 
DEFINE A2 = CRPT ACCE 
CRPT ACCE = SUM(A1,A2) 

♦GOTO 200 
♦LABEL 70 
$ 

$ REACTIONS 
$ 

*IF("NM" NE REAC) : ♦GOTO 90 
♦GOTO 200 
♦LABEL 90 
$ 

$ STRESSES 
$ 

♦IF("NM" NE STRE): ♦GOTO 110 
$ FORM [S] [DQ/DV] 
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★DCALL (TR , CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VSTRE CNAME=STRE Q=QX 

★ XQT DCU 

CHANGE 1 CRPT STRE 1 1 CRPT STR1 1 1 
$ FORM [DS/DV] [Q] 

★ XQT ES 
RESET OPER-T 

IDQ = "XLIB" SEL STRESS 
U = "XLIB'* VIBR MODES 1 1 1 , "NMODES" 

T = TMAT VSTRE 
★XQT AUS 

DEFINE SO = "XLIB" TMAT VSTRESS 
DEFINE SI = TMAT VSTRESS 

TMAT DSDV = SUMO'OVDX" SI "MOVD" SO) 

★DC ALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSDV CNAME=STRE Q=QX 

★ XQT DCU 

CHANGE 1 CRPT STRE 1 1 CRPT STR2 1 1 
$ FORM S [D PHI / DV] [Q] 

★XQT ES 
RESET OPER=T 
IDQ = "XLIB" SEL STRESS 
U * DXDV AUS " IDV" 1 1 , "NMODES" 

T = TMAT DSTRE 

★DC ALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSTRE CNAME=STRE Q=QX 
★XQT AUS 

DEFINE STR1 = CRPT STR1 

DEFINE STR2 = CRPT STR2 

DEFINE STR3 = CRPT STRE 

TMP = SUM(STR1 , STR2) 

CRPT STRE = SUM(TMP , STR3) 

★GOTO 200 
★LABEL 110 
★GOTO 300 
★LABEL 200 
★XQT DCU 

CHANGE 1 CRPT "NM" 1 1 DXDV "NM" "IDV" 1 
COPY 1 "XLIB" DXDV "NM" "IDV" 1 
PRINT "XLIB" DXDV "NM" "IDV" 1 
★LABEL 300 
! IIB = IIB + 1 
★ JGZ , -1 (NNB , 10) 

! IIB = FREE ( ) : ! NNB = FREE ( ) : ! LIB1 = FREE ( ) : ! LIB2 = FREEO 
! LIB3 = FREEO : ! TNAME = FREEO : ! CNAME = FREEO 
! Q = FREEO 
★END 

Runstream TR DBACK 4 

Runstream TR DBACK 4 recovers sensitivities in the fixed-mode semianalytical method with 
the mode acceleration method. 

s 

$ (TR DBACK 4) - BACK TRANSFORMATION FOR DERIVATIVES USING 
$ SEMIANALYTICAL METHOD 

$ WITH THE MODE ACCELERATION METHOD 

$ 

★XQT AUS 

DEFINE X = "XLIB" VIBR MODE 1 1 1, "NMODES" 

DEFINE E = "XLIB" VIBR EVAL 1 1 
DEFINE DKX = DKX AUS 
DEFINE DMX = DMX AUS 
DEFINE ROMG = "XLIB" ROMG AUS 

DEFINE XTDX = "XLIB" XTDX 

DEFINE XOMD = "XLIB" XOMD 

DEFINE XOME = "XLIB" XOME 

★IF("IDSP" NE 0): TABLE(NI="NMODES" , NJ="NMODES") : RDDX 
$ CALCULATE VELOCITY TERM 


X0M1 = CBR(XOME.RDKX) 

XOMK = CBD(XDM1 ,R0MG) 

! NJDM = TOC,NJ("XLIB" XTDX MASK MASK MASK) 

! NBDM * TOC, NINJ( "XLIB" XTDX MASK MASK MASK) 

♦IF("NJDM" NE "NBDM") : XOKC = CBR (XOMK , XTDX) 

♦IF("NJDM" EQ "NBDM"): XOKC = CBD (XOMK , XTDX) 

XOMC = CBR(XOME.RDDX) 

XQD = SUM (XOKC , -1.0 XOMC) 

$ CALCULATE ACCELERATION TERM 
DKXO = CBD(DKX , RDMG) 

APPL FORC 887 = SUM(DKXO, -1.0 DMX) 

$ CALCULATE DERIVATIVE OF THE PSEUDOSTATIC TERM 
DEFINE USTAT = "XLIB" STAT DISP 1 1 
FSL1 = PROD (DKDV, USTAT) 

APPL FORC 888 = SUM (DFDV , -1 . 0 FSL1) 

♦XQT SSOL 

RESET SET=887 , KLIB="XLIB" , KILIB="XLIB" , REAC=0 
♦XQT SSOL 

RESET SET=888 , KLIB="XLIB" , KILIB="XLIB" , REAC=0 

$ 

$ LOOP OVER ALL RESPONSE QUANTITY TYPES 

$ 

! IIB * 1 i ! NNB = NBCK 
♦LABEL 10 
♦XQT AUS 

! NM * DS, 11 IIB" , 1,1 ("XLIB" BACK LIST 1 1) 

! I ERR = TOC,IERR("XLIB" SEL "NM" MASK MASK) 

♦JNZ(IERR,300) 

$ 

$ DISPLACEMENTS 
$ 

♦IF ( "NM" NE DISP): ♦GOTO 30 
DEFINE IDJK = "XLIB" SEL DISP 
DEFINE DAC1 = STAT DISP 887 
DEFINE DUST = STAT DISP 888 
DEFINE XTDX = "XLIB" XTDX 
DEFINE DACC = "XLIB" TMAT DACC 
TMAT DUST = SVTRAN( IDJK , DUST) 

TMAT DAC1 = SVTRAN ( IDJK ,DAC1 ) 

$ VELOCITY TERMS 
TMAT DVL1 = SVTRAN (IDJK, XQD) 

$ 

♦DCALL (TR , CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DUST CNAME=DISP Q=A 
T0CC(1 CRPT DISP 1 1) : N2=DSP1 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DVL1 CNAME^DISP Q=QX1 
TOCCCl CRPT DISP 11): N2=DSP2 

♦DCALL (TR, CRPT) LIB1*"XLIB" LIB2="XLIB" LIB3=1 TNAME=DVEL CNAME*DISP Q^QXl 
TDCCd CRPT DISP 1 1) : N2*DSP3 

♦DCALL (TR, CRPT) LIB1«"XLIB" LIB2=1 LIB3="XLIB" TNAME=DAC1 CNAME^DISP Q=QX2 
TDCCd CRPT DISP 11): N2=DSP4 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=*DACC CNAME=DISP Q=QX2 
TDCCd CRPT DISP 11): N2=DSP5 
♦XQT U4 
VU 

DEFINE D1 = CRPT DSP1 
DEFINE D2 = CRPT DSP2 
DEFINE D3 = CRPT DSPS 
DEFINE D4 = CRPT DSP4 
DEFINE DS = CRPT DSPS 

CRPT DISP = SUM (D1 , D2 , -1.0*D3, D4, -l.O^DS) 

♦GOTO 200 
♦LABEL 30 
$ 

$ VELOCITIES 
$ 

♦IF("NM" NE VELO): ♦GOTO 50 
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♦DCALL (TR , CRPT) LIB1="XLIB" LIB2="XLIB" LIB3-1 TNAME=VVEL CNAME-VELO Q=QX1 
♦GOTO 200 
♦LABEL 50 
$ 

$ ACCELERATIONS 

$ 

♦IF( "NM" NE ACCE): ♦GOTO 70 

♦DCALL (TR , CRPT) LIB1-"XLIB" LIB2="XLIB" LIB3=1 TNAME=VACC CNAME=ACCE Q=QX2 
♦GOTO 200 
♦LABEL 70 
$ 

$ REACTIONS 

$ 

♦ IFC'NM" NE RE AC) : ♦GOTO 90 
♦GOTO 200 

♦LABEL 90 

$ 

$ STRESSES 

$ 

♦ IFC'NM" NE STRE): ♦GOTO 110 

$ 

$ FORM [S] [DU/DV] 

$ 

♦LIBS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦XQT ES 
RESET OPER-T 
IDQ = SEL STRESS 

U= M XLIB" STAT DISP 888 : T = "XLIB" TMAT DTM1 

U="XLIB" XQD AUS : T * "XLIB" TMAT DTM2 

U="XLIB" STAT DISP 887 : T = "XLIB" TMAT DTM4 

♦LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
♦XQT AUS 

♦DCALL (TR, CRPT) LIB1-"XLIB" LIB2=1 LIB3="XLIB" TNAME=DTM1 CNAME-STRE Q=A 
T0CCC1 CRPT STRE 11): N2=STR1 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3-"XLIB" TNAME=DTM2 CNAME-STRE Q=qXl 
T0CC(1 CRPT STRE 11): N2=STR2 

♦DCALL ( TR , CRPT ) LIB1="XLIB" LIB2-"XLIB" LIB3=1 TNAME-SD CNAME-STRE Q=qXl 
T0CC(1 CRPT STRE 11): N2=STR3 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3= H XLIB" TNAME=DTM4 CNAME-STRE q=qX2 
T0CC(1 CRPT STRE 1 1) : N2-STR4 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME-SP CNAME=STRE Q=qX2 
T0CC(1 CRPT STRE 11): N2-STR5 

$ 

$ FORM [DS/DV] [U] 

$ 

♦XQT ES 
RESET OPER-T 
IDQ = "XLIB" SEL STRESS 
U = "XLIB" STAT DISP 11 : T = TMAT SF 

U = "XLIB" XOMD AUS : T = TMAT SD 

U = "XLIB" XOME AUS : T = TMAT SP 

♦XQT AUS 

DEFINE SO = "XLIB" TMAT SF : DEFINE SI = TMAT SF 
TMAT DSF = SUMO'OVDX" SI "MOVD" SO) 

DEFINE SO * "XLIB" TMAT SD : DEFINE SI = TMAT SD 

TMAT DSD = SUMO'OVDX" SI "MOVD" SO) 

DEFINE SO = "XLIB" TMAT SP : DEFINE SI = TMAT SP 

TMAT DSP = SUMO'OVDX" SI "MOVD" SO) 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSF CNAME=STRE Q=A 
T0CC(1 CRPT STRE) : N2=STR6 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSD CNAME=STRE Q=QX1 
T0CC(1 CRPT STRE) : N2-STR7 

♦DCALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSP CNAME-STRE Q=QX2 
T0CC(1 CRPT STRE) : N2-STR8 
♦XQT U4 
VU 
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DEFINE SI = CRPT STR1 
DEFINE S2 = CRPT STR2 
DEFINE S3 = CRPT STR3 
DEFINE S4 * CRPT STR4 
DEFINE S5 = CRPT STR5 
DEFINE S6 = CRPT STR6 
DEFINE S7 = CRPT STR7 
DEFINE S8 = CRPT STR8 

CRPT STRE = SUM (SI , S2, -1.0*S3, S4, -1.0*S5, S6 , -1.0*S7, -1.0*S8) 

♦GOTO 200 
♦LABEL 110 
♦GOTO 300 
♦LABEL 200 
♦XQT DCU 

CHANGE 1 CRPT "NM" 1 1 DXDV "NM" "IDV" 1 
COPY 1 "XLIB" DXDV "NM" "IDV" 1 
PRINT "XLIB" DXDV "NM" "IDV" 1 
♦LABEL 300 
! IIB = IIB + 1 
♦ JGZ, -1 (NNB , 10) 

! IIB = FREE ( ) : ! NNB = FREEC) : ! LIB1 = FREE() : ! LIB2 = FREE ( ) 

! LIBS = FREE ( ) : » TNAME = FREE() : ! CNAME = FREEC) 

! Q = FREEC) 

♦END 

Runstream TR DPHI 3 

Runstream TR DPHI 3 implements the modified modal method for calculating eigenvector 
derivatives and is called from sensitivity calculation runstream TR DXDV 6. 

$ 

$ (TR,DPHI , 3) - CALCULATE EIGENVECTOR DERIVATIVES USING THE 
$ MODIFIED MODAL METHOD 

$ 

♦XQT AUS 

INLIB=10 : OUTLIB-IO 
DEFINE MO - "XLIB" "MNAME" 

DEFINE DK = 1 DKDV SPAR 
DEFINE DM * 1 DMDV SPAR 

DEFINE XO = "XLIB" VIBR MODE 1 1 1, "NMODES" 

DEFINE AJK = AJK 

DEFINE EV = "XLIB 11 VIBR EVAL 1 1 
MX = PRDD(DM,X0) 

AKK = XTYDC-.5 XO,MX) 

TABLE CNI*"NMODES" , NJ="NMODES" ) : AJK 

TABLECNI-l , NJ*"NMODES") : UNIT : J*1 , "NMODES" : 1.0 

! J = 1 : ! NJ * NMODES 
♦LABEL 10 

! EJ = DS , " J" , 1 , 1 ("XLIB" VIBR EVAL 1 1) 

! MEJ = -EJ 

DEFINE XJ = "XLIB" VIBR MODE 1 1 "J","J" 

DKDM = SUM (DK, "MEJ" DM) 

MX = PROD (DKDM, XJ) 

DLAM = XTY (X J ,MX) 

! DLAM - DS, 1,1,1(10 DLAM AUS 1 1) 

AF1 = PROD ("DLAM" MO, XJ) 

AF2 = SUMCAF1 -1.0 MX) 

♦IF(" J" EQ 1): 11 APPL FORC = UNI0N(AF2) 

♦IF(" J" NE 1): 11 APPL FORC = UNION, U(AF2) 

AA = XTY(XO,MX) 

DENI * SUMC-1.0 EV, "EJ" UNIT) 

DEN 2 = PROD (EV, DENI) 

TABLE, U : DEN 2 : I="J" : J=1 : 1.0 
FACT = RECIPCDEN2) 

AAB = PROD("EJ" FACT.AA) 

DEI : OPER=XSUM : DEST,U=AJK AUS 
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SOURCE=AAB : JS=1 : JD= M J" : EX1 

SQURCE=AKK : IS="J" : JS=1 : ID="J" : JD="J" : EX1 
! J = J + 1 

♦ JGZ , -1 (NJ , 10) 

DXDV = CBR(X0 , AJK) 

♦XQT SSOL 

RESET KLIB="XLIB" , KILIB="XLIB" , QLIB=11, REAC=0, EP=0 

♦ XQT AUS 

DEFINE D1 * 10 DXDV AUS 
DEFINE D2 = 11 STAT DISP 
DXDV AUS "IDV" = SUM (D1 , D2) 

* DELETE 10 
* DELETE 11 
♦RETURN 
♦END 


Runstream TR CRPT 

TR CRPT is a utility runstream which performs the transformation from modal to physical 
basis for a single response quantity at a set of critical times. It is heavily used by the TR DBACK 
runstreams. 

$ 

$ (TR, CRPT) - FORM CRITICAL POINT RESPONSE TABLE 

$ 

♦0NLINE=0 
♦XQT AUS 


NCRIT = TOC , NI ( "LIBl" CRTI 

"CNAME" 

1 

1) 

$ 

NUMBER 

OF 

TIMES 

ND * T0C,NI("LIB2" TMAT 

"TNAME" 

1 

1) 

$ 

NUMBER 

OF 

RESP. QUANTITIES 

NQ = TOC, NJ ( "LIB2" TMAT 

"TNAME" 

1 

1) 

$ 

NUMBER 

OF 

MODES 


! NJQ = TOC , NJ ("LIB3" "Q" AUS MASK MASK) 

! ISTP = 0 

$ LOOP OVER ALL RESPONSE QUANTITIES 
! I = 1 
! N = ND 

INLIB = 21 : OUTLIB = 21 
♦LABEL 20 
DEI 

S0URCE="LIB2" TMAT "TNAME" 

ID = 1 : IS* "I" 

DEST-TDNE "TNAME" "I" 1 
EX1 

TABLE(NI»"NQ" , NJ="NCRIT" ) : XBAR CRIT "I" 

$ LOOP OVER NUMBER OF CRITICAL POINTS 
! II = 1 
! NN = NCRIT 
♦LABEL 40 
DEI 

! TIME = DS , "II" , 1 , "I" ("LIB1" CRTI "CNAME" 1 1) 

! ISTP = TIME/DT + .5 
! ISTP = ISTP + 1 
! IB = ISTP/NJQ 
! 1ST = IB^NJQ 

♦ IFO'IST" NE "ISTP"): ! IB = IB + 1 
! J * IB - 1 * NJQ : ! J = ISTP - J 

SOURCE * "LIB3" "Q" AUS MASK MASK "IB" : JS="J" : DEST,U=XBAR CRIT "I" 
JD * "II" : EX1 
! II = II + 1 
* JGZ , -1 (NN ,40) 

DEFINE T = TONE "TNAME" "I" 

DEFINE XBAR = XBAR CRIT "I" 

Cl AUS "I" = RPROD(T.XBAR) 

TOCC(CI AUS "I") : NJ=1 
DEFINE Cl = Cl AUS "I" 

♦IF ("I" EQ 1): 1 CRPT "CNAME" = UNION(CI) 

♦ I F ( " I " GT 1): 1 CRPT "CNAME" = UNION, U(CI) 
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! 1 = 1 + 1 
♦JGZ,-1(N,20) 

♦0NLINE=1 

! NCRIT = FREEO : ! ND = FREE() : ! NQ = FREEO ! NJq “ FREEO 

! ISTP * FREEO : ! IB = FREEO : ! 1ST = FREEO 

♦END 

Runstream TR DIAG 

TR DIAG is a utility runstream that solves the reduced-order eigenproblem based on a given 
set of basis vectors to uncouple a reduced system. 

$ 

$(TR DIAG) - MAKE AN ARBITRARY VECTOR SET M AND K ORTHONORMAL 

$ 

♦ XQT AUS 

0UTLIB=10 : INLIB=10 

DEFINE K*1 K: DEFI M=1 "MNAME" 

DEFINE X = "VLIB" VIBR MODE 
I JCODE=10000 

KX=PRQD(K ,X) : SYN K 10000 "NMODE" = XTYS(X.KX) 

MX-PR0D(M , X) : SYN M 10000 "NMODE" = XTYS(X,MX) 

! ZER0=NM0DE-1 

♦ JZ (ZERO ,1003) 

♦XQT STRP 

RESET SOURCE- 10 , DEST=10 

♦ JGZ (ZERO, 1004) 

♦ LABEL 1003 

♦ XQT AUS 

0UTLIB=10: INLIB=10 

!K-DS 2 1 1(10 SYN K MASK MASK) 

!M“DS 2 1 1(10 SYN M MASK MASK) 

!EVAL=K/M 

TABLE (NI-1 , NJ=1 ) : SYSEVEC: J-l: 1.0 
TABLE(NI®1 ,NJ=1) : SYS EVAL : J*l: "EVAL 1 ' 

♦ LABEL 1004 

♦ XQT AUS 

0UTLIB=10: INLIB=10 
DEFINE E=SYS EVEC 
DEFI X = "VLIB" VIBR MODE 
X ORTH 1 1=CBR(X,E) 

DEFINE X=X ORTH 1 1 

"VLIB" VIBR MODE 1 1 “UNION (X) 

DEFINE E=SYS EVAL 

"VLIB" VIBR EVAL 1 1=UNI0N(E) 

♦END 
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